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INTRODUCTION 

In  the  past  2i  decades  the  approach  outlined  by  St.  Denis  and 
Pierson  (Reference  1)*  for  the  estimation  of  oscillatory  ship  motion 
statistics  in  irregular  seas  has  become  firmly  established  in  engineer¬ 
ing  practice.  This  approach  applies  strictly  only  to  ship  responses 
which  can  be  assumed  to  be  linear  functions  of  wave  height.  The 
credibility  of  the  approach  was  established  by  exercising  experimen¬ 
tally  as  many  as  possible  of  the  mathematical  consequences  of  the 
basic  mathematical  model.  In  particular  it  was  found  possible  to 
a)  synthesize  observed  responses  in  irregular  waves  in  both  time  and 
frequency  domains  by  means  of  responses  obtained  in  regular  waves, 
and  b)  identify  the  fundamental  linear  response  functions  from  obser¬ 
vations  in  both  irregular  and  transient  waves  (Reference  2)*. 

In  general,  when  non-linear  responses  become  of  importance 
there  is  no  agreed  universal  model  for  dealing  with  the  irregular 
sea  case.  However,  when  the  non-linearities  may  be  considered 
"weak"  in  some  sense  one  of  the  conceptual  approaches  which  have 
been  proposed  has  considerable  attraction.  This  is  the  functional 
series  model.  Among  the  attractions  are  that  the  model  is  suitable 
for  any  reasonably  well  behaved  wave  input  (regular,  transient  or 
random)  and  since  the  model  contains  the  completely  linear  system 
as  a  special  case  it  appears  to  be  a  logical  extension  of  present 
practice.  In  addition,  prediction  methods  for  scalar  response 
spectra  are  available  and  it  appears  that  the  statistics  of  maxima 
may  be  approximated.  Finally,  it  is  possible  to  closely  relate  the 
functions  required  by  the  model  to  deterministic  hydromechanical 
analyses  and  experiment  because  the  effects  of  hydrodynamic  "memory" 
which  complicate  the  usual  analysis  are  automatically  accounted  for. 


1.  St.  Denis,  M,  and  Pierson,  W.J.,  "On  the  Motions  of  Ships  in 

Confused  Seas,"  SNAME  Vol.  61,  1953. 

«Hf 

2.  Dalzell,  J.F.,  "The  Input-Output  Approach  to  Seakeeping  Problems: 

Review  and  Prospects,"  T.SR.  Symposium  $-3,  Seakeeping  1 953~ 1 973 » 
Society  of  Naval  Architects  and  Marine  Engineers,  October  1973. 
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References  3*  through  7*  represent  a  contracted  bibliography  of  some 
of  the  fundamental  mathematical  background  to  the  model  and  References 
8*  and  9*  discuss  some  of  the  potential  for  application  to  seakeeping 
problems. 

In  practical  use  the  functional  series  must  be  truncated  at 
some  degree  of  non-linearity,  thus  becoming  a  functional  polynomial. 
Within  the  last  decade  the  functional  polynomial  model  of  degree  two 
has  been  applied  with  some  success  to  the  non-linear  ship  responses 
in  the  un-restored  modes  of  motion,  in  particular  to  added  resistance, 
slow  drift  oscillations,  etc.  By  drawing  upon  the  analytical  background 


3.  Weiner,  N.,  "Non-Linear  Problems  in  Random  Theory,"  The  Technology 

Press  of  MIT  and  John  Wiley  and  Sons,  Inc.,  1958. 

& 

4.  Barrett,  J.F.,  "The  Use  of  Functionals  in  the  Analysis  of  Non- 

Linear  Physical  Systems,"  Journal  of  Electronics  and  Control, 

15,  No.  6,  December  1963. 

•ft 

5.  George,  O.A.,  "Continuous  Non-Linear  Systems,"  Doctoral  Disser¬ 

tation,  Department  of  Electrical  Engineering,  MIT,  July  1959. 

6.  Ku,  Y.H.,  and  Wolf,  A. A.,  "Vol terra-Wei ner"  Functionals  for  the 

Analysis  of  Non-Linear  Systems,"  Journal  of  Franklin  Insitutute, 
281 ,  No.  1,  January  1966. 

'fc 

7.  Bedrosian,  E.  and  Rice,  S.O.,  "The  Output  Properties  of  Vol terra 

Systems  (Non-Linear  Systems  with  Memory)  Driven  by  Harmonic 
and  Gaussian  Inputs,"  Proceedings  of  the  IEEE,  Vol.  59,  No.  12, 
December  1971. 

"fc 

8.  Vassi lopoulos,  L.A.,  "The  Application  of  Statistical  Theory  of 

Non-Linear  Systems  to  Ship  Motion  Performance  in  Random  Seas," 
Ship  Control  Systems  Symposium,  Annapolis,  November  1966. 

* 

9.  Bishop,  R.E.D.,  Burcher,  R.K. ,  and  Price,  W.G.,  "The  Uses  of 

Functional  Analysis  in  Ship  Dynamics,"  Proceeding,  Royal 
Society  of  London,  A.  332,  1973. 
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afforded  by  References  10*  through  13*  it  was  possible  to  develop 
practical  means  of  doing  cross-bi-spectral  analyses  of  towing  tank 
observations  of  added  resistance,  Reference  14*,  and  following  this 
a  demonstration.  Reference  15*,  that  a)  the  non-linear  added  resis¬ 
tance  frequency  response  function  can  be  identified  from  irregular 
wave  experiments  via  cross-bi-spectral  analysis,  b)  the  mean  added 
resistance  and  the  spectrum  of  resistance  can  be  synthesized  from 
l  the  linear  and  non-linear  frequency  response  functions,  c)  that 

time  histories  of  added  resistance  can  be  synthesized  according  to 
the  functional  polynomial  model,  and,  d)  that  hydrodynamic  theory 
can  be  developed  for  the  required  frequency  response  functions 
I  (Reference  16)*. 


10.  Tick,  L.J.,  "The  Estimation  of  the  "Transfer  Functions"  of 

Quadratic  Systems,"  TECHNOMETRICS,  3,  No.  4,  1961. 

* 

11.  Hasselman,  K. ,  "On  Non-Linear  Ship  Motions  in  Irregular  Waves," 

JSR  J_0 ,  No.  1,  1966. 

12.  Shaman,  Paul,  "Bi-Spectral  Analysis  of  Stationary  Time  Series," 

Scientific  Paper  #18,  Statistical  Laboratory,  School  of 
Engineering  and  Sciences,  N.Y.U.,  January  1964. 

*13.  Rosenblatt,  M.  and  Van  Ness,  J.W.,  "Estimates  of  the  Bi-Spectrum 
of  Stationary  Random  Processes,"  Technical  Report  11,  Nonr 
562 (29) /1 1 ,  Division  of  Applied  Mathematics,  Brown  University, 
Providence,  R.I.,  March  1964. 

14.  Dalzell,  J.F.,  "Cross-Bi-Spectral  Analysis:  Application  to 

Ship  Resistance  in  Waves,"  Journal  of  Ship  Research,  Vol.  18, 
No.  1,  March  1974,  pp.  62-72. 

*15.  Dalzell,  J.F.,  "Application  of  the  Functional  Polynomial  Model 
to  the  Ship  Added  Resistance  Problem,"  Eleventh  Symposium  on 
Naval  Hydrodynamics,  University  College,  London,  1976. 

16.  Dalzell,  J.F.  and  Kim,  C.H.,  "An  Analysis  of  the  Quadratic 
Frequency  Response  for  Added  Resistance,"  Journal  of  Ship 
Research,  Vol.  23,  No.  3,  September  1979. 


I 


I 


3 


TR-2275 


It  is  becoming  clear  that  even  when  the  functional  polynomial 
model  is  not  explicitly  assumed  the  results  of  many  investigations  of 
the  second  order  forces  on  ships  (References  17*  and  18*  for  example) 
are  compatible  with  the  functional  polynomial  model  of  the  second 
degree.  It  is  equally  clear  that  this  "1 i near+quadrat i c"  model  will 
not  suffice  for  non-linear  response  problems  where  the  non-linearity 
affects  the  amplitude  at  excitation  frequency  of  oscillatory  response 
to  harmonic  excitation.  (In  the  “linear  plus  quadratic"  or  second 
degree  system  the  non-linearities  only  produce  new  response  frequency 
components.)  Thus  it  is  logical  to  consider  applying  the  functional 
polynomial  of  degree  three  to  some  seakeeping  problems.  It  is 
realistic  to  consider  whether  this  is  necessary  since  the  theoretical 
hydrodynamici st  would  consider  the  non-linearities  of  the  added  (cubic) 
degree  to  be  “third-order".  Despite  the  implications  of  the  words 
“third-order",  “third-order"  forces  on  heaving  cylinders  have  been 
found  (Reference  19)*  not  to  be  negligibly  small  for  large  heave 
amplitudes  and  in  fact  magnitudes  approaching  those  of  "first-order" 
or  linear  forces  have  been  reported.  The  work  of  References  20*  and  21* 


17.  Newman,  J.N.,  "Second  Order  Slowly  Varying  Forces  on  Vessels  in 
Irregular  Waves,"  International  Symposium  on  the  Dynamics 
of  Marine  Vehicles  and  Structures  in  Waves,  University  College, 
London,  April  1974. 

JL 

’l8.  Pinkster,  J.A.  and  van  Oortmerssen,  G.,  "Computation  of  the  First 
and  Second  Order  Wave  Forces  on  Bodies  Oscillating  in  Regular 
Waves,"  Proceedings  of  the  Second  International  Conference 
on  Numerical  Ship  Hydrodynamics,  University  of  California, 
pp.  136-156,  September  1977. 

19.  Tasai,  F.  and  Koteratama,  W.,  "Non-Linear  Hydrodynamic  Forces 

Acting  on  Cylinders  Heaving  on  the  Surface  of  a  Fluid,"  Reports 
of  the  Research  Institute  for  Applied  Mechanics,  Vol.  XXIV, 

No.  77,  1976. 

£ 

20.  Dalzell,  J.F.,  "A  Note  on  the  Form  of  Ship  Roll  Damping,"  SIT-DL- 

76-1887,  Davidson  Laboratory,  Stevens  Institute  of  Technology, 
May  1976,  (Also:  Journal  of  Ship  Research,  Vol.  22,  No.  3, 
September  1978). 

* 

21.  Dalzell,  J.F.,  "Estimation  of  the  Spectrum  of  Non-Linear  Ship 

Rolling:  The  Functional  Series  Approach,"  S I T-0L- 76- 1 89^ ,  Dav¬ 
idson  Laboratory,  Stevens  Institute  of  Technology,  May  1976, 
AD-A031  055/7G1. 
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has  shown  that  if  some  of  the  long  accepted  ideas  about  non-linear 
ship  rolling  are  to  be  accommodated  by  the  functional  polynomial 
approach,  the  polynomial  must  be  at  least  of  degree  three.  Thus 
there  was  encouragement  in  prior  work  to  consider  further  the  applic¬ 
ability  of  the  third  degree  model  to  seakeeping  problems. 

The  main  avenue  of  establishing  credibility  in  the  cases  of 
linear  and  1 i near+quadrat i c  systems  was  through  the  capability  of 
decomposing  observations  in  irregular  seas  by  analysis.  In  the 
linear  case  this  has  been  through  cross-spectral  analysis.  In  the 
1 i near+quadrat i c  case  this  has  been  through  cross-bi-spectral 
analysis.  In  the  1 i near+quadrat ic+cubic  case  (the  functional  poly¬ 
nomial  of  degree  three)  the  analysis  approach  is  not  established. 

The  objectives  of  the  present  work  were  to  further  explore 
the  applicability  of  the  third  degree  functional  polynomial  to 
seakeeping  problems  and  to  attempt  the  development  of  an  analysis 
approach  by  which  third  degree  non-linearities  in  the  response  of 
ships  to  random  waves  might  be  interpreted  and  characterized. 
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THE  FUNCTIONAL  SERIES  MODEL 


Genera  1 i t i es 

The  mathematical  applications  involved  in  the  present  work  were 
first  expressed  by  Wiener,  Reference  3,  and  have  since  been  slowly 
amplified  and  simplified  for  engineering  purposes  to  a  point  that 
specialized  text  book  coverage  has  begun  to  appear  (Reference  22*  for 
example).  All  the  fundamental  mathematical  development  has  taken 
place  in  the  context  of  system  or  communication  theory  where  typically 
a  "system"  has  a  single  "input"  or  excitation  and  a  single  "output" 
or  response.  Clearly,  the  general  seakeeping  problem  is  more  com¬ 
plicated  since  in  that  case  the  physical  ship  system  is  conceptually 
excited  by  elemental  wave  components  arriving  simultaneously  from  a 
continuous  range  of  headings'-the  ship  in  short  crested  seas  may  at 
the  very  least  be  considered  a  multiple  input  system.  For  present 
purposes  however  the  single  excitation  assumption  must  be  made,  and 
this  translates  into  the  case  of  long  crested  seas.  In  the  event  that 
the  ship  has  forward  speed,  additional  specializations  must  be  made 
in  order  that  the  available  mathematical  background  be  utilized.  These 
are  that  the  mapping  of  wave  frequency  into  encounter  frequency  be 
essentially  single  valued,  and  that  the  excitation  be  defined  at  a 
point  which  translates  with  the  ship.  These  specializations  correspond 
to  what  is  usually  done  in  towing  tank  experiments  as  well  as  the 
specializations  customary  in  hydromechanical  analyses.  Essentially, 
the  restrictions  imposed  by  the  available  systems  background  suggest 
applicability  to  a  useful  (but  not  all  inclusive)  sub-class  of  sea¬ 
keeping  problems. 

The  Functional  Series 

Within  the  class  of  problems  so  defined  it  will  be  possible  to 
identify  a  response,  Y(t),  of  interest  (a  function  of  time,  t)  and 
an  excitation,  X(t),  which  is  typically  identified  as  a  wave 


* 

22.  Rugh,  W.J.,  "Nonlinear  System  Theory;  The  Volterra/Wiener  Approach," 
Johns  Hopkins  University  Press,  1981. 
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elevation.  Thus  for  present  purposes  it  is  reasonable  to  assume  that 
the  excitation,  X(t),  is  a  zero  mean  function  whether  it  be  random  or 
deterministic.  Thus,  following  References  14,  15,  21,  as  a  first  step 
it  is  hypothesized  that  the  response,  Y (t),  is  a  sufficiently  regular 
function  that  it  may  at  least  be  expanded  in  an  infinite  functional 
ser i es : 


» 


Y(t)  =  I 
n*o 


r 

L9 


„(t  ,t 

n  i  ; 


t  )X  (t-t, ) • • *X( t-t  )dt  dt  •••dt  i  (1 

n  l  n  i  2  1i 


(Omission  of  limits  on  integrals  here  and  throughout  this  report  signify 
limits  of  -»  and  +°°. ) 


Each  term  in  the  series  as  written  above  is  homogeneous  func¬ 
tional  of  degree  n.  The  terms  are  said  to  be  homogeneous  because  a 
change  in  X ( t )  to  c  X(t)  (where  c  is  constant)  results  in  multiplication 

of  the  term  of  degree  n  by  (cn).  The  kernels,  g  (t  •••t  ),  are  "time 

n  i  n 

^  invariant",  since  they  are  considered  to  be  functions  only  of  time 

d? fferences .  In  Che  present  application,  the  wave  system  varies  with 
time,  but  not  the  ship.  Consequently  the  properties  of  the  ship  are 
contained  wholly  within  the  kernels. 

)  Reference  4  indiciates  that  the  series  expansion  is  unique  if 

all  kernels  are  completely  symmetrical  in  the  variables;  that  is 


I 


9n(V 


g„(t  ,t 
n  2 


for  any  rearrangment  of  the  variables 
results  from  this  restriction.) 


3 

t . 
J 


)... 


(No  loss  of  generality 


(2) 


According  to  Reference  6,  the  functional  series  converges  for 
bounded  excitation  so  long  as  the  sum  of  integrals  of  the  absolute 
value  of  all  kernels  is  less  than  +■>.  The  same  is  true  with  stochastic 
excitation,  if  in  addition  the  input  is  strictly  stationary  with 
bounded  moments  of  all  orders.  One  physical  way  of  looking  at  the 
convergence  restrictions  on  the  kernels  is  that  the  "memory"  of 
the  ship  must  not  be  enormously  long,  a  restriction  intuitively 
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acceptable  for  most  if  not  all  seakeeping  responses.  The  restrictions 
on  the  excitation  are  those  which  have  been  tacitly  accepted  for  some 
time  in  seakeeping  research. 

So  long  as  the  series,  Equation  1,  converges,  its  value  is  that 
of  the  kernel  of  zeroth  degree  (gQ)  when  the  excitation,  X(t),  is  zero. 
In  the  present  context,  responses  in  the  absence  of  waves  are  not  of 
interest  so  that  the  first  term  in  the  series  may  be  dropped  in  all 
subsequent  development. 

Impulse  and  Frequency  Response  Functions 

The  kernels  in  Equation  1  may  be  considered  as  describing  the 
system  through  a  series  of  nC^  degree  impulse  response  functions.  It 
is  presumed  that  each  impulse  response  function  is  sufficiently  smooth 
and  integrable  so  that  there  is  no  trouble  about  existence  of  an  n-fold 
Fourier  transform.  Accordingly,  it  is  assumed  that  to  each  n1^1  degree 
impulse  response  function  there  corresponds  an  nth  degree  frequency 
response  function,  Gn (wj ,u>2* *un) .  The  transform  pairs  relating  impulse 
and  frequency  response  functions  may  be  defined  as  follows: 


gn(t  ' 

n  x  2 


- 


<2n)n 


G  (w  ,w 
n  !  2 


•  ) 

n 


Exp 


L  J-1 


,.n. 


u> .  t  -  doi  dcu  « 
J  Jj  1  2 


(3) 


G  (uj  ,  w  •  •  *0)  ) 

n  i  2  n 


g  (t  ,t  »»*t  )  Exp  -i  I  w.t.  dt  dt 
n  1  2  n  L  ja|  J  U  12 


(wh 


ere  the  to. 

J 


are  ci rcular 


frequencies) . 


•dt 

n 


The  first  degree  frequency  response  function  is  the  familiar 
linear  one.  The  second  degree  frequency  response  function  is  the  one 
treated  in  References  1A,  15,  and  16  in  an  application  to  added  resis¬ 
tance.  Regardless  of  the  degree,  the  basic  importance  of  the  trans¬ 
form  of  the  impulse  response  function  is  the  same;  that  is,  convolution 
in  the  time  domain  usually  corresponds  to  multiplication  in  the  fre¬ 
quency  domain,  and  (perhaps  more  importantly  in  the  present  context) 
the  typical  hydromechanical  seakeeping  analysis  ends  up  in  the 
frequency  domain. 
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As  a  consequence  of  the  assumed  symmetry  of  the  impulse  response 
functions  and  the  transform,  Equation  4,  the  nC^  degree  frequency 
response  function  is  also  symmetric  in  its  arguments.  That  is, 

G  (w  ,ui  •  •  •  uj  )  *  G  ,ui  •••u  )••• 
n  i  2  n  n21  n 

for  any  and  all  rearrangements  of  the  u>. .  Additionally,  because  the 
impulse  response  functions  are  real: 


G  (-^i  ,—j»  ,■ 
n  1  2 


G  (ui  ,»■  •••«  ) 
n  1  2  n 


where  the  star  denotes  the  complex  conjugate,  and  all  arguments  on  the 
left  hand  side  are  negative. 
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THE  THIRD  DEGREE  FUNCTIONAL  POLYNOMIAL 

It  is  clear  from  the  form  of  Equation  1  that  the  complexity  of  any 
answer  which  might  result  would  increase  geometrically  with  the  order 
of  the  term.  Consequently,  it  is  hoped  that  the  non-linearities  in 
any  given  system  are  weak  enough  that  the  series  may  be  truncated  at 
a  relatively  few  terms,  in  which  case  it  is  termed  a  "functional  poly¬ 
nomial".  The  objective  of  the  present  work  is  to  investigate  the 
applicability  of  the  functional  polynomial  of  the  third  degree  to 
seakeeping  problems.  Truncating  Equation  1  at  n  *  3  results  in  the 
fundamental  time  domain  mathematical  model  for  present  purposes: 


g  (t  ) X( t-t  ) dt 
11  11 


+  g  (t  ,t  )X(t-t  )X( t-t  ) dt  dt 

J  J  2  1  2  1  2  12 

+  I  [  [  g  (t  ,t  ,t  )X(t-t  ) X ( t-t  )X( t-t  )dt  dt  dt 

J  J  J  3  1  2  3  1  2  3123  (5) 

This  model  is  the  same  as  that  of  References  14  and  15  with  the  func¬ 
tional  of  third  degree  added.  Though  it  may  seem  physically  and 
intuitively  justifiable  to  accept  convergence  of  the  series,  Equation  1, 
the  acceptance  of  the  functional  polynomial,  Equation  5>  as  an  adequate 
engineering  «,p  yroximat  ion  hinges  purely  upon  how  well  it  works  in 
practice  with  whatever  problem  is  at  hand.  The  series  itself.  Equation  1, 
is  thought  applicable  only  to  weakly  nonlinear  systems.  If  the  system 
responds  critically  to  a  given  level  of  excitation  (as  for  instance 
a  "jump"  or  an  instability  in  the  sense  of  the  Hill  equation)  many 
terms  in  the  series  would  be  necessary  for  even  a  rough  approximation 
and  practical  application  is  doubtful.  On  the  other  hand  (apart  from 
the  added  res i stance /s  low  dr  if t  phenomina  which  are  quadratic  non- 
linearities)  inmany  seakeeping  situations  where  nonlinearities  might 
be  expected  a  purely  linear  treatment  often  yields  results  which 
appear  within  reason,  and  the  model  of  Equation  5  may  represent  an 
improvement. 
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DETERMINISTIC  IDENTIFICATION  OF  FREQUENCY  RESPONSE  FUNCTIONS 

Given  the  polynomial.  Equation  5,  and  the  Fourier  transform 
definition.  Equation  A,  there  are  three  frequency  response  functions 
of  interest  in  the  present  problem: 

(Linear)  ■*  Gj  (w) 

(Quadratic)  -*■  G  (u,,ai 2) 

(Cub i c)  -*•  G  2  (uj  j  ,  o>2  t  ^3  ) 

The  interpretation  of  G^  (u>)  is  of  course  identical  to  that  of  pure 
linear  theory,  it  is  simply  the  normalized  amplitude  and  phase  of 
linear  steady  state  system  response  to  sinusoidal  excitation.  The 
quadratic  response,  G2 (w,  ,u)2) ,  has  been  interpreted  in  References  14  and 
15,  and  what  nominally  remains  is  to  interpret  the  cubic  frequency 
response  function.  In  words,  G2  (o>j  ,±w2)  expresses  the  normalized 
steady  state  response  at  frequencies  (utj+u^)  to  two  cosinusoids  of 
frequencies  w  and  u>2,  due  to  interaction  of  the  two  frequency  com¬ 
ponents.  The  units  of  G2(<Dj,w2)  are  (response  uni  t)/(exci  tation  unit)2. 
Similarly,  G3 (u  ,±u>2 ,±w3)  expresses  the  normalized  steady  state  response 
at  frequencies  (u^tu  ±a>3)  to  three  cosinusoids  of  frequencies  o>  ,  w2 
and  w  ,  due  to  interaction  of  the  three  frequency  components.  The 
units  of  G3  (ui1  ,w2,io3)  are  (response  uni  t)/(exci  tation  unit)3. 

When  the  polynomial  of  third  degree,  Equation  5,  is  considered 
to  be  the  system  there  arise  some  serious  complications  not  present 
in  the  second  degree  model  considered  in  References  14  and  15.  For 
this  reason  it  is  thought  best  to  clarify  the  contents  of  the  last 
paragraph  in  stages  which  correspond  to  "experiments"  of  increasing 
complexity. 

If  the  system  is  completely  linear  [g2(t i«t2)  and  93(tj,t2,t3) 

*  0  in  Equation  f]  experimental  characterization  takes  the  form  of 
exciting  the  system  with 

X ( t )  ■  A  cos(wt  -e)  (6a) 

and  interpreting  the  steady  state  response  as: 

Y(t)  -  A  Re[G,(a>)  Exp  (itui  -  iefl  (6b) 
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from  which  the  linear  frequency  response  function  G^w)  may  be  derived 
by  a  frequency  analysis.  It  should  be  noted  that  e  is  an  arbitrary 
phase  angle  inserted  in  Equation  6a  to  reflect  the  situation  which  is 
normally  inherent  in  the  reduction  of  observations  from  a  physical 
experiment.  In  linear  theoretical  analyses  much  the  same  process  is 
followed  except  that  "A"  is  customarily  taken  to  be  unity,  "e"  as 
zero,  and  the  excitation  written  as 

X  ( t )  *  Exp  [jtu] 

with  the  understanding  that  the  real  part  of  the  resulting  complex 
solution  for  Y(t)  is  to  be  taken.  The  result  is  identical  in  form 
to  Equation  6b,  and  the  theoretical  expression  for  G^u)  may  be 
extracted.  It  should  be  remarked  before  proceeding  further  that 
the  customary  theoretical  representation  of  excitation  as  a  complex 
quantity  with  the  real  part  of  the  resulting  response  "understood" 
is  a  completely  valid  procedure  only  for  linear  systems.  Application 
of  this  procedure  to  nonlinear  systems  of  the  present  type  results 
in  an  incomplete  answer  at  best.  In  order  to  produce  theoretical 
analyses  compatible  with  physical  experiment  for  the  present  nonlinear 
system  the  excitation  must  be  carried  as  an  explicitly  real  quantity, 
or  as  sums  of  complex  quantities  and  their  conjugates. 

Returning  to  the  nonlinear  system.  Equation  5,  the  simplest 
"experiment"  is  to  excite  the  system  with  a  single  cosinusoid  just  as 
in  the  characterization  of  the  fully  linear  system.  Thus,  as  before: 

A 

X ( t)  ■  Excitation  *  A  cos  (wt  -  e)  (7a) 

Expressing  Equation  7  a  as  (A/2)  [[Exp  ( i  tio  -  e)  +  Exp  (-itu>  +  eQ, 
substitution  into  Equation  5,  appl i cation  of  the  Fourier  transform 
definition.  Equation  A,  and  consideration  of  symmetries,  results  in 
the  following  expression  for  the  response  o?  the  nonlinear  system  to 
the  single  cosinusoid.  Equation  7a: 

Response  ■  Y(t)  ■  Re 

+  j  A2Re^2 (w,-u>)J 

+  j  A2RejG,(uj,w)  Exp(it2w-  i 2e )^J  (7b) 

+  jf  A3RefG3  («,w,w)  Exp(it3w 


^  A3G3  (u>,w,-w)  }Exp(i  tw  -  ie)J 


12 


TR'2275 


The  compl ication  mentioned  earlier  is  immediately  apparent  in  the  first 
term  of  Equation  7b.  The  cubic  nonlinearity  may  produce  response  at 
the  excitation  frequency  as  well  as  at  the  third  harmonic  of  the  exci¬ 
tation.  While  the  values  of  the  functions  G0(u>,u)  and  G2(w,-w),  and 
the  values  of  G,(w,u,u)  may  in  principle  be  inferred  by  a  frequency 
analysis  of  a  single  experiment,  more  than  one  experiment  is  required 
to  separate  G  j  (uj)  and  G3  (w,u),-<d)  .  The  complication,  however,  may  not 
arise  in  theoretical  analysis  because  G3(w,uj,-lo)  is  the  component 
,  which  varies  as  excitation  amplitude  cubed. 

Continuing  the  discussion  of  Equation  7b,  the  quantity  most 
reduced  and  extracted  in  seakeeping  experiments  in  waves  corresponds 
to  the  ratio  of  amplitude  of  response  at  frequency  u>  to  amplitude  of 

>  a  regular  wave  of  frequency  ui.  In  terms  of  Equation  7b  this  ratio  is: 

R(ui)  =  |G^  (u)  +  A2G3  (w,w,-w)  |  (8) 

When  Gg  (w.w.-u)  is  zero,R(w)  is  [  G  J  (u»)  |  and  is  invariant  with  excitation 

>  amplitude  as  would  be  expected  for  a  system  of  second  degree.  When 

the  cubic  nonlinearity  is  present  the  amplitude  can  vary  with  excitation 
amplitude  in  a  variety  of  ways,  since  in  general  both  G^u)  and  G3(w,u),-u>) 
are  complex  and  this  must  be  taken  account  of  prior  to  performing  the 

>  absolute  value.  If  the  ratio  of  real  and  imaginary  parts  of  G1(u)  is 

the  same  as  the  corresponding  ratio  for  G3  (u.to.-w) ,  R (uj)  ,  (an  "equivalent" 
linear  response  ratio)  could  be  quadratic  in  excitation  amplitude. 

Because  the  single  cosinusoid  experiment  involves  only  specialized 
^  portions  of  the  quadratic  and  cubic  frequency  response  functions  it  is 

clear  that  more  complicated  "experiments"  must  be  carried  out  to  define 
the  remaining  portions,  and  these  take  the  form  of  dual  and  triple 
"tone"  experiments  where  the  excitation  is  taken  as  the  sum  of  2  and 
^  3  cosinusoids. 

Because  of  the  cubic  term  the  dual  tone  experiment  outlined  in 
Reference  15  takes  a  more  complicated  form.  In  particular  the  excitation, 
I  X ( t ) ,  is  assumed  to  be  composed  of  two  cosinusoids  of  different  amplitudes 

Aj  and  A2,  and  different  frequencies  and  u>2  as  shown  at  the  top  of 


» 
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Table  1.  In  order  to  be  consistent  with  physical  experiments,  arbitrary 
phases  Ej  and  e2  are  included.  The  response,  Y(t),  to  the  dual  tone 
excitation,  X ( t )  is  shown  in  the  lower  part  of  Table  1.  The  result  is 
derived  exactly  as  indicated  for  Equation  7b,  the  assumed  excitation 
is  substituted  in  Equation  5  in  complex  form,  the  transform  definition, 
Equation  4,  is  applied,  and  as  much  advantage  as  possible  is  taken  of 
symmetry.  A  further  complication  about  the  response  at  the  excitation 
frequencies  is  evident  from  the  first  three  terms  in  the  expression  for 
Y(t),  Table  1.  To  illustrate,  the  response  at  frequency  Uj  is  written 
out  as  fol lows : 

Re[{A1G1  (uj)  +  -^  A^G3  (u)l  ,oj1  ,-Uj^)  +  j  A^A2G3  (w2 , -u)2,u1 )  }Exp(  i  tw1  -  iEjf] 

The  first  two  terms  inside  the  curly  brackets  are  the  same  as  in  the 
corresponding  term  of  the  single  cosinusoid  experiment,  the  third 
represents  a  third  degree  interaction  between  the  two  frequency  com¬ 
ponents.  The  fourth  through  seventh  lines  of  the  response  in  Table  1 
are  the  same  as  the  quadratic  terms  described  in  Reference  15- 
The  eighth  line  involves  the  third  harmonics  of  the  two  excitation 
frequencies  as  would  be  expected  from  the  single  tone  experiments. 
Equation  7b.  The  last  four  lines  involve  a  frequency  of  a  new  form; 
that  is, 

2u.  ±  a), 
i  J 

and  the  corresponding  special  values  of  quadratic  frequency  response 
function,  G 3 (co .  ,oj.  ,±o>j) . 

Table  2  indicates  the  response  of  the  system  of  Equation  5  to 
a  three  tone  excitation.  In  this  case  all  the  types  of  frequency 
components  noted  in  the  two  tone  experiment  reappear.  Essentially, 
the  results  of  Table  1  are  repeated  for  all  the  possible  combinations 
of  two  of  three  frequencies.  The  last  four  lines  of  the  response 
shown  in  Table  2  involve  four  new  frequencies: 

(oJj+U^+UJj) 

(h,l+u,2*W3) 

(uj  -w  +u>  ) 

1  2  3 

<WV 
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TABLE  1 


RESULTS  OF  THE  TWO  COSINUSOID  EXPERIMENT 


Excitation  *  X ( t )  -  Ajcos  (uijt  -  +  A, cos  (w:t  -  ej 

Response  ■  V(t) 

2  "1 

«  l  RejjAj.Gi  (uj)  +  T  Aj63(wj*“j*"u’j)}  ExP('twj  ”  iej)J 

+  2  AlA2ReR Exp ( i tw,  -  ie2)J 

+  \  A2AxRe[^3  (‘4i2»“CLi2  ,U‘l)  Exp^itU1l  ‘  iel|] 

+  lAlG2(“l-*U1)  +  T  A2G2U2'"W2) 

♦  l  y  AjRejs^Uj.w.)  Exp  (ft2u>j  -  f2cj)J 

+  AlA2Re[G2(w1,(i)2)  Exp  {it(u>1  +  «2)  -  i  (e1  +  e2)}J 
+  A1A2Re[c2(u.2f-w1)  Exp  {it(u»2  -  u^)  -  i(e2  -  e^J 

+  ^  ^  AjRe^  j  (uij  ,Wj  »<*)j  )  Exp  (it3wj  *  i3€j)^ 

+  |>  AfA2R«j63(u1(ultu2)  Exp  {it(2wl  +  w,)  -  it(2Cj  +  e2)}J 
+  f  AiA2Re[® 3  (t*>i  ,0*! , -w2)  Exp  { i  t(2uij  -  u^)  -  it(2ej  -  ej  )J 

+  i  AiA2Re[f  ^  (u»2  ,a»2  ,U>1  ^  ExP  * '  +  "  't(2e2  +  e1)}J 

+  i  AiA|Re[^  3  (u»2  ,u»2 ,  -o)j )  Exp  { i  t  (2u»2  -  u>t)  -  it(2e,  - 
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Exci tat  ion 
Response 


TABLE  2 

RESULTS  OF  THE  THREE  C0SINUSI00  EXPERIMENT 

:  3 

X(t)  ■  I  A. cos  (tui.  -  C .  ) 


j-1 


J  J 


Y(t) 
3 


I  R«|{A.G  (w.)  ♦  t  A?G3(u>.',u>.,-ui.)  }  Exp(itiu.  -  ic.) 
j-lLJiJ  J  J  J 

2  J,  Jk  Exp(lt»k  -  ujj 


j l,  T*js2<»j--‘jl 


I  T  AfRefG. («.,«.)  Exp( i t2w.  -  i2e.f| 
j-1  *  J  L‘  J  J  J  JJ 

jll  AJAk**Eal“k»“j)  ^^“k  +  “j>  *  }<*k-e.)}] 

j-!  AjV*[®2(“k'*“j)  ‘“j5  •  it(tk  ’  cj)}] 

3  r  i 

X  AjR«|G3(uj,Uj,uij)  Exp(it3uij  -  I3*j)J 

I  I  \  A?A.  Reffi,  («.,«.  ,w.)  Exp{  t  t(2u.  ♦  w.)  -  f  (2c .  ♦  c.  )  fl 
lc-1  ji*k  *  J  *  L3  J  J  k  J  *  J  k  J 

3  3 


*  I  I  r  A?A,  Re|G,(w.,tii.,-u>v )  Exp{it(2u>.  ■  w.  )  •  i(2c.  -  e.)}l 

k-1  j#k  J  *  L 3  J  J  *  J  *  J  * 

+  \  AlA2A3R*[®3^“l,l*'2,“3^  *  <*>2  *  <*>3 )  >Exp{- i  ( tj  *  tj  ♦  Cj)}^ 

^  E*p{it(w1  *  uij  -  uij)  }Exp{- i  (tj  ♦  *2  •  c3)H 

^  A1A2A3A*  ^3  ’  ”*“2 '*"3  ^  "  m2  ^  “3)  )Exp{-i  (Cj  *  E,  *  E  3 )  >  I 

*  4  A  A  A  ReIG  (-ui  )  Exp{  i  t(w  ♦  ui.  -  *  )  }Exp{-i  (c  ♦  c,  •  c  ) ) ! 

*  1*3^3  14.  3  3fc  I  3  i 
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The  first  is  the  triple  sum  frequency  which  is  associated  with  the  value 
of  G3  (ojj  ,u>2  ,(d3)  in  the  positive  octant  of  tr  i -frequency  space  (and  by 
symmetry  the  negative  octant).  The  three  remaining  frequencies  are 
associated  with  the  cubic  frequency  response  functions  in  the  remaining 
six  octants  of  tr i -frequency  space. 

Thus  in  principle  the  three  tone  experiment,  repeated  for  enough 
combinations  of  the  three  frequencies,  can  completely  define  the  cubic 
frequency  response  function.  If  the  frequencies  chosen  are  incommen¬ 
surate  a  very  precise  frequency  analysis  can  in  principle  allow  the 
identification  of  G,  (u^  ,w2  ,ui3 )  by  picking  out  only  the  response  at 
the  sum  frequency.  This  is  the  attitude  expressed  in  the  systems 
literature.  Reference  22  for  example. 

The  results  in  Table  2  suggest  how  an  analytical  identification 
of  the  various  frequency  response  functions  may  be  made.  Essentially 
the  analytical  problem  amounts  to  the  consideration  of  the  response 
to  three  superimposed  cosinusoids,  and  the  comparison  of  the  resulting 
complex  coefficients  of  each  of  the  various  time  factors  with  the 
corresponding  terms  in  Table  2. 

In  the  experimental  case  of  interest  here  (towing  tank  experi¬ 
ments)  finite  tank  length  tends  to  limit  the  overall  length  of  experi¬ 
ment,  and  the  incommensurate  frequency  approach  might  well  be  very 
difficult.  If  it  is  possible  to  choose  excitation  frequencies  such 
that  a  complex  periodic  wave  form  results,  relatively  simple  fre¬ 
quency  analysis  might  serve  as  a  data  reduction  procedure.  To  indicate 
the  dimension  of  the  problem  in  a  more  concrete  way,  all  the  expected 
frequencies  in  the  response  to  1 ,  2  and  3  cosinusoids  have  been  listed 
in  Table  3.  In  the  single  tone  experiment  there  are  U  expected  fre¬ 
quencies,  in  the  two  tone  case  there  are  13,  and  in  the  three  tone 
case  32. 

Tables  1  through  3  together  suggest  how  the  deterministic 
identification  problem  might  be  approached.  The  first  step  would 
be  to  run  single  tone  experiments  to  isolate  Gj (w) ,  G2(w,-w),  G2(w,w), 

G3  (us,u>,-u>) ,  and  G3(u),oj,w)  as  noted  in  the  discussion  of  Equation  7b 
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jjk 


TABLE  3 

FREQUENCIES  PRESENT  IN  THE  RESPONSE  TO  THE  p  COSINUSOID  EXPERIMENT 


TR-2275 


(Experiments  with  more  than  one  excitation  amplitude  would  be  required 
to  separate  G^Cu)  and  G  (u,w, -w) . )  In  principle,  once  these  special 
cases  are  identified  it  would  be  possible  to  progress  to  the  two  tone 
experiment.  It  will  be  noted  from  Table  3  that  in  the  two  tone  exper¬ 
iment  seven  of  the  13  frequency  components  involve  values  of  response 
functions  derivable  from  the  single  tone  experiment.  Thus  the  concen¬ 
tration  would  be  upon  the  quadratic  sum  and  difference  frequencies 
and  the  cubic  terms  arising  from  two  tone  interactions.  The  object 
of  the  two  tone  experiment  would  thus  be  to  identify  G  oj2), 

G2(iolt-w2),  the  special  values  of  the  cubic  response  function  of  the 

form  G.  (10.  ,u>.  ,+w.  ) ,  and  to  separate  the  excitation  frequency  response 
j  j  J  k 

G3  (toj ,  -OK  ,0)^)  from  the  linear  and  cubic  response  obtained  in  the  single 
tone  experiment.  Again  in  principle,  the  one  and  two  tone  experiments 
would  serve  to  identify  all  but  the  last  four  frequency  components  of 
the  three  tone  experiment  so  that  there  would  be  a  reduction  of  the 
problem  of  selecting  three  excitation  frequencies  such  that  all  the 
32  frequencies  shown  for  the  three  tone  experiments  would  be  distinct. 
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FREQUENCY  DOMAIN  SIMULATION 
OF  A  CUBIC  SYSTEM 

In  order  to  proceed  much  further  with  the  present  work  it  was 
desirable  to  produce  a  specific  example  of  a  system  of  the  third 
degree  in  order  to  see  what  qualitative  features  a  cubic  nonlinearity 
might  introduce  into  the  response  to  random  waves,  and  to  enable  the 
simulation  of  the  "data"  which  would  be  necessary  in  a  later  attempt 
at  identification.  A  reasonable  course  of  action  was  suggested  by 
the  work  of  Reference  21  where  a  modified  form  of  the  classic  single 
degree  of  freedom  roll  equation  was  expanded  in  a  functional  series. 

Thus  for  present  purposes  it  was  assumed  that  the  response, 

Y ( t) ,  and  the  excitation,  X(t),  of  the  simulated  system  were  related 
by  the  following  differential  equation: 

3  r—  .  ~ 

l  A  (Y(t)}J  +  B  .  {Y  ( t)  }J  +  C  -  {Y  ( t)  }J  =  X(t)  (9) 

j  =  1  LJ  J  J  J 

where  A^ ,  Bj  and  Cj  are  constants  and  A^,  Bj  and  C:  are  not  zero  so  that 
the  system  can  contain  a  significant  linear  response.  In  Reference  21 
the  constants  corresponding  to  j=2  were  taken  to  be  zero  so  that  no 
quadratic  nonlinearities  were  present,  the  equation  was  expanded  in  a 
functional  series  to  fifth  degree  and  the  frequency  response  functions 
were  derived  by  the  i ncommensurate  frequency  technique  of  Reference  7. 

For  present  purposes  an  expansion  up  to  the  third  degree  was  wanted 
and  thus  most  of  the  work  required  for  the  present  case  is  documented 
in  Reference  21.  Applying  the  incommensurate  frequency  approach  of 
Reference  7  to  Equation  9,  and  truncating  the  results  at  the  third 
degree  to  be  consistent  with  the  cubic  model,  Equation  5,  results  in 
the  following  expressions  for  the  simulated  linear,  quadratic  and  cubic 
frequency  response  functions: 

GjU)  =  1/D1  ( iu>)  (10) 

G^u^.u^)  *  -D2(-w1w2)  01(w2)  GjUj  +  u,2)  (11) 
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G  j  («-, )  Gj  (u.^ ) 


(12) 

where  the  auxiliary  function  is  defined: 

D  (a)  =  A  a2  +  B  a  +  C  03) 

n  n  n  n 

Some  remarks  are  perhaps  in  order  before  proceeding  further.  It 
may  be  noted  in  Equations  11  and  12  that  the  expressions  for  the  quad¬ 
ratic  and  cubic  frequency  response  functions  involve  ail  the  response 
and  auxiliary  functions  of  lesser  degree.  Unfortunately,  this  continues 
on  for  frequency  response  functions  of  higher  degree  so  that,  despite 
the  fact  that  the  exponent,  j,  of  Equation  9  is  limited  to  3,  frequency 
response  functions  of  all  degrees  may  be  derived  for  Equation  9.  In 
effect.  Equations  10  through  12  are  not  a  complete  solution.  This  is 
not  too  bothersome  in  the  present  context  since  the  postulated  model, 
Equation  5,  assumes  that  all  response  functions  of  greater  degree  than 
three  are  zero,  and  since  all  that  is  wanted  here  are  a  set  of  response 
functions  which  have  the  required  symmetry,  and  which  may  be  evaluated. 

Given  the  coefficients  of  Equation  9,  numerical  evaluation  of 
Equations  10  through  12  is  straightforward.  What  was  wanted  for  the 
present  work  was  one  fairly  reasonable  system,  and  the  selection  of  the 
coefficients  was  partly  arbitrary  and  partly  trial  and  error.  The 
coefficients  Aj,  B3,  were  arbitrarily  fixed  as  follows: 

A>  =  1/(2tt)2 

Bj  *  1/4tt 

Ci  -  1 

These  choices  define  a  simple  linear  single  degree  of  freedom  system  in 
which  the  response  amplitude  ratio  is  unity  at  zero  frequency  and  2.0 


G3  (uij  ,U)„  ,0J3)  =  “G^UJj  +  U)?  +  ^,)  J-Dj  (-  iuljU.'-.uJj)  G 

+  y  D0(-uy(uj„  +  oJ 3 ))  G,(w2>u)3)  G}  (uy) 

+  y  +  ^3))  G„  (uy  ,uy  )  G  1  (co2  ) 

+  \  D  (-to  (uj  +  GJ„))  G„  (u>.  .a,.)  G  (us  ) 


3 


3  '  1 


1  '~3' 
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at  resonance.  The  resonant  frequency  is  2^  radians/second.  Thus  the 
linear  subsystem  was  made  generally  similar  to  a  lightly  damped  heave 
response  to  waves  or  a  relatively  heavily  damped  roll  response  to  wave 
slope.  The  implied  time  scale  (resonant  period  of  1  second)  is  typical 
of  small  model  testing.  It  may  be  noted  in  Equations  10  through  12 
that  if  is  zero,  and  G,  .,w  )  will  also  be  zero 

for  any  choice  of  and  u2.  for  the  numerical  work  to  be  later 
described  it  was  thought  to  be  of  advantage  to  make  sure  that  the 
quadratic  and  cubic  frequency  response  functions  would  go  to  zero  if 
the  absolute  value  of  any  frequency  argument  exceeded  some  cutoff  value. 
Accordingly,  in  the  numerical  evaluation  of  Gj (oj) ,  provision  was  made 
to  multiply  the  computed  value  of  G  (u)  in  the  range  7^  <  |u|  <  8tt  by  a 

*  A 

factor  decreasing  linearly  from  unity  to  zero.  In  effect,  G  (w)  was 
computed  in  accordance  with  Equation  10  in  the  range  -711  <  w  <  7t,  the 
computed  value  outside  this  range  was  attenuated  to  zero  at  | a> |  =  8ir, 
and  Gj  (to)  was  set  equal  to  zero  for  J u> |  >  8tt. 

The  next  step  was  to  select  the  quadratic  coefficients  A2.  B0,  C2, 
in  Equation  9.  For  lack  of  better  information,  the  intention  w.y$  !j 
produce  a  quadratic  frequency  response  function  qualitatively  similar 
to  those  shown  in  References  15  and  16.  Qualitatively,  the  functions 
in  those  references  are  double  humped  in  the  bi -frequency  plane.  Gener¬ 
ally,  the  maximum  absolute  value  of  G2(w,w)  defines  the  maximum  of  one 
hump  and  G2(w,-w)  the  maximum  of  the  other.  The  function  tends  to 
zero  for  bi-frequency  (0,0)  and  to  be  small  for  large  values  of  either 
frequency  component.  Given  the  values  of  A  ,  B  ,  C1  established  pre- 
viously,  G2(w,w)  and  G2(w,-gj)  were  evaluated  in  accordance  with  Equ¬ 
ation  11  for  trial  and  error  selections  of  A2,  B„  and  C2.  The  first 
conclusion  was  that  C2  should  be  zero  if  G2(0,0)  was  to  be  zero.  With 
an  initial  selection  of  A2  and  B2  which  yielded  reasonable  qualitative 
behavior  of  G2(w,w)  and  G2(oj,-w),  the  entire  function  was  computed.  It 
was  found  that  the  A^  coefficient  produced  relatively  significant  values 
of  response  at  large  values  of  frequency  argument  (|wj|  J4  |w2|)  which 
made  the  function  qualitatively  dissimilar  to  that  desired,  and  this 
coefficient  was  finally  also  zeroed.  The  final  selection  of  quadratic 
coefficients  was  as  follows: 
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A;  =  C.  =  0 
B:  =  -0.265  A: 

This  choice  yielded  the  overall  qualitative  behavior  desired,  a  maximum 
value  of  G ,  (oj ,  — 0, )  of  about  uni  ty  at  -j  =  2ti,  and  a  maximum  value  of 
|G,  (ui.w)  |  of  -0.33  at  the  same  frequency.  Noting  Equation  7b,  for 
single  tone  excitation  of  unity  ampl i tude,  these  values  imply  a  shift 
in  the  mean  of  about  1/2  and  a  second  harmonic  of  about  1/6. 

Given  the  linear  and  quadratic  frequency  response  functions,  the 

coefficients  A  ,  B  ,  C  were  to  be  selected  in  order  to  complete  the 

3  3  3 

simulated  system.  In  this  case  there  was  no  previous  qualitative 
guidance  available  since  no  mapping  of  the  cubic  frequency  response 
function  for  any  system  was  available.  It  seemed  reasonable  to  select 
these  coefficients  upon  the  basis  of  the  single  tone  response.  Equation  7b 
that  is,  to  compute  j  G 1  (  2tt)  +  0.75  G3  (2tt,2tt,-2tt)  |  ,  which  is  the  ampli¬ 
tude  of  the  response  at  the  resonant  frequency  to  unity  amplitude 
excitation  at  resonance.  The  objective  was  to  select  coefficients  which 
would  make  this  result  about  A/3,  (2/3  of  the  linear  response  amplitude) 
on  the  assumption  that  this  might  yield  an  "equivalent'1  linear  frequency 
response  resembling  something  which  might  have  been  experienced  in  ship 
motions  experiments.  Initial  trial  and  error  numerical  work  showed 
that  the  A3  coefficient  had  to  nearly  vanish  in  order  to  achieve  the 
desired  result.  It  was  suspected,  as  in  the  quadratic  case,  that  a 
non  zero  A3  coefficient  would  produce  significant  values  of  G  (coj , u>2, ou^) 
for  high  frequency  arguments.  A  preliminary  mapping  of  the  function 
with  non-zero  A3  tended  to  confirm  this  and  the  final  values  of  coeffi¬ 
cients  selected  were: 


B3  =  0.0032  Bj 

c  3  =  0.0075  C3 

With  this  selection  the  maximum  absolute  value  of  G3(w,u),-w)  was  0.92 
and  the  maximum  absolute  value  of  G3(uj,uj,w)  was  about  0.11;  that  is, 
the  third  harmonic  response  to  single  tone  excitation  would  be  about 
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1/9  the  amplitude  of  the  contribution  to  the  fundamental  frequency. 

This  last  is  not  inconsistent  with  ship  motions  experiments  where 
significant  third  harmonic  response  is  rare  or  un-noticed.  The  non¬ 
zero  C3  coefficient  produces  a  slight  static  softening  effect  which 
is  also  not  inconsistent  with  ship  hydrostatics. 

As  has  been  previously  implied,  one  of  the  main  motives  for  the 
present  work  has  been  the  experimental  situation  where  the  normalized 
response  amplitude  is  not  invariant  with  excitation  amplitude  as  it 
should  be  for  a  purely  linear  system.  It  was  thus  of  interest  to 
evaluate  the  normalized  response  ("equivalent"  linear  response)  for 
the  simulated  system.  R(oj)  of  Equation  8  was  evaluated  for  a  number 
of  excitation  amplitudes  ("A"  of  Equation  7a)  and  the  results  are 
plotted  in  Figures  1  and  2.  The  evaluations  were  carried  out  at  fre¬ 
quency  steps  of  V10,  the  actual  results  are  plotted  as  points  and 
straight  lines  are  drawn  between.  The  amplitude  range  chosen  was 
A  *  0  to  5.  The  somewhat  surprising  range  of  the  results  suggested 
the  presentation  in  Figures  1  and  2.  The  frame  at  the  top  of 
Figure  1  covers  the  range  0  ^  A  <  2,  the  frame  at  the  bottom  extends 
the  range  to  3,  and  Figure  2  extends  the  amplitude  range  to  4  and  5. 

It  is  clear  from  the  figures  that  the  cubic  nonlinearities  which  would 
be  noticeable  in  an  experiment  with  the  simulated  system  are  concen¬ 
trated  about  the  linear  resonance.  Normalized  response  above  w  =  12 
is  essentially  invariant  with  excitation  amplitude.  That  below 
w  :  3.5  is  weakly  dependent  on  amplitude.  The  decrease  of  the  normal¬ 
ized  response  near  resonance  is  about  as  might  be  expected  from  the 
method  of  choosing  coefficients,  up  to  an  excitation  amplitude  of  2. 

Above  this  amplitude  the  trend  of  response  reverses  and  the  normalized 
response  grows  steadily  with  excitation  amplitude.  The  radical  trends 
shown  in  Figure  2  probably  have  no  qualitative  parallel  in  known  sea¬ 
keeping  experiments.  On  the  otherhand  the  qualitative  and  quantitative 
behavior  shown  at  the  top  of  Figure  1  for  an  excitation  amplitude  range 
up  to  about  1.0  has  been  seen  quite  a  lot,  and  it  is  conceivable 
physically  that  nonlinearities  could  produce  inflections  in  the  normalized 
response  at  a  given  frequency,  as  implied  in  the  lower  frame  of  Figure  1. 
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Two  interim  conclusions  arise.  First,  the  simulated  system 
may  be  reasonably  representative  of  some  nonlinear  seakeeping  situations 
so  long  as  the  magnitude  of  excitation  amplitude  is  two  or  less.  Second, 
the  cubic  model  seems  capable  of  reflecting  far  more  complicated  quali¬ 
tative  behavior  than  had  been  imagined  at  the  outset. 
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SIMULATED  IMPULSE  RESPONSES 

The  next  step  in  the  simulation  was  to  invert  the  frequency 
domain  simulation  into  the  time  domain;  that  is,  to  derive  digital 
linear,  quadratic  and  cubic  impulse  response  functions.  The  time 
domain  operation  intended  was  to  carry  out  a  discrete  version  of 
Equation  5  which  may  be  written  as  follows: 

Y(t)  -  Y(n)  -  Y1(n)  +  Y2(n)  +  Y3(n)  (WO 

where: 

Yl(n)  *  l  g  j  X ( n  -  j)  (15) 

j  J 

Y2(n)  =  l  l  g]  X(n  -  j)  X(n  -  k)  (16) 

j  k  J 


Y3(n)  *  l  l  l  9?,  X(n  -  j)  X(n  -  k)  X(n  -  l)  (17) 

j  k  t  JKi 

where  in  terms  of  Equation  5: 

n  *  t/At 
j  -  tx/At 
k  =  t2/ At 
l  -  tj/At 

At  ■  a  sampling  interval 

Y(n)  *  Y(nAt)  ■  the  sampled  response 
X(n)  *  X(nAt)  *  the  sampled  excitation 

and  the  digital  impulse  response  functions  are: 


9j  a  9  !  (j A t) 

•  At 

(18) 

9jk  *  92(j4t' 

kAt)  •  At2 

(19) 

9jk*  ■  93(j4t- 

kAt,  J.At)  •  At2 

(20) 

(for  integer  j,  k  and  l) . 
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The  domain  of  the  summations  in  Equations  15  through  17  is  over 
non-zero  values  of  the  digital  impulse  response  functions.  The  sampling 
interval  must  be  chosen  small  enough  that  the  implied  trapezoidal  inte¬ 
grations  are  reasonable,  and  for  the  system  defined  in  the  previous 
section  a  sampling  interval  of  0.0625  *  At  was  chosen. 

Computation  of  the  impulse  response  functions,  as  in  Reference  15, 
amounts  to  carrying  out  numerically  the  applicable  Fourier  transform, 
(Equation  3)  for  each  of  the  frequency  response  functions  for  a  range 
of  j,  k  and  l  sufficient  to  allow  eventual  truncation  of  the  function. 

The  required  numerical  operations  were  carried  out  with  the  Fast 
Fourier  Transform.  In  particular  let: 

2n  p 
“l  "  At  N 

2  TT  £ 
w2  *  At  N 

2  it  r 
U3  *  At  N 

and  assume  trapezoidal  integration  of  Equation  3,  with  integer  values  of 
p,  q  and  r  within  the  following  ranges: 

-N/2  <  p  <  N/2 
-N/2  <  q  <  N/2 
-N/2  <  r  N/2 

1  l  Gj  (2irp/ AtN)  Exp(i2irjp/N)  (21) 

P 

~  I  I  G  (2:rp/AtN,  2irq/AtN)  Exp{  i2ir  ( jp+kq)/N}  (22) 

N2  P  q  2 

~  1  I  I  G,(2TTp/AtN,2irq/AtN,2Trr/AtN)  Exp{  i 2it (jp+kq+Hr)/N)  (23) 
N3  p  q  r 

In  these  equations  numerical  values  of  G^w),  G2(ultu>2)  and  G3  (Wj  ,w2  ,u>3) 
were  computed  in  accordance  with  Equations  10  through  13  and  the  values 
of  the  coefficients  of  Equation  9  which  were  given  in  a  previous  section. 


Then: 


3jk  = 


’  j  k£ 
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Equation  21  is  exactly  what  the  inverse  Fast  Fourier  Transform 
does,  apart  from  the  leading  factor  of  1/N,  and  that  values  of  Gj (— lu) 
must  appear  as  if  aliased  in  the  second  half  of  the  complex  N  point  FFT 
frequency  domain  array.  An  additional  numerical  complication  is  that 
the  inverse  FFT  yields  a  "circular"  time  function  which  corresponds 
properly  to  positive  and  negative  values  of  time  only  if  the  frequency 
function  presented  to  it  is  zero  for  H/k  <  jp|  <  N/2.  With  the  provision 
discussed  previously  that  G1  (w)  =  0  for  |m|  >  8k,  and  the  value  of  At 
selected,  this  requirement  is  satisfied  for  N  *  128,  which  was  the  value 
used  throughout  the  computations.  Once  the  inverse  FFT  is  performed 
gj  is  defined  for  (-N/2  +  1)  <  j  <  N/2.  The  value  for  j  =  -N/2  is  not 
defined,  but  the  omission  is  not  crucial  since  if  the  computation  para¬ 
meters  have  been  selected  properly  the  function  must  be  negligible  at 
the  ends  of  the  j  range. 

The  computation  of  Equation  22  is  equivalent  to  first  doing 
Equation  21  N  times  for  constant  q,  {(-N/2  +  l)  <  q  <  N/2}.  This  yields 
a  partial  transform  of  )  which  is  a  function  of  j  and  q 

(tx  and  u  ),  and  N  additional  inverse  transforms  with  respect  to  q(o>2) 
complete  the  computation. 

The  computation  of  Equation  23  is  equivalent  to  doing  Equation  22 
N  times  for  the  range  of  r,  which  yields  a  mixed  function  of  j ,  k  and  r 
(t  ,  t2,  and  w3),  and  N2  inverse  transforms  with  respect  to  r(w3)  complete 
the  work.  It  should  be  noted  in  this  connection  that  with  N  =  128  the 
computer  memory  requirements  for  Equations  21  and  22  are  quite  modest. 

To  accomplish  Equation  23  with  all  elements  in  memory  at  once  requires 
of  the  order  of  four  million  words  of  memory  so  that  in  the  present 
instance  the  computation  of  Equation  23  was  made  in  stages  with  storage 
of  intermediate  results. 

As  a  practical  matter  minimization  of  the  necessary  range  of 
j,  k  and  l  is  advantageous  in  the  application  of  the  discrete  impulse 
responses,  Equations  15  through  17.  The  results  of  the  computations  of 
Equations  21,  22  and  23  were,  in  all  three  cases,  quite  small  relative 
to  the  maximum  for  negative  values  of  j,  k  and  l.  This  is  the  result 
expected  from  the  form  of  Equation  9.  Since  the  coefficients  are 
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fixed  constants,  the  system  should  be  realizable  in  real  time.  This 

implies  that  g,(t  ),  g„(t  ,t  )  and  g  (t  ,t  ,t  )  should  be  exactly  zero 

11  -  :  ; 

if  any  time  argument  is  negative,  because  negative  time  arguments  imply 
that  the  "future"  must  be  known  in  order  to  predict  the  "present".  The 
fact  that  the  values  of  the  computed  impulse  responses  were  negligible 
but  not  exactly  zero  is  probably  attributable  to  the  truncation  of 
G^gj)  described  previously  and  to  rounding  error  in  the  computations. 

Truncation  of  the  computed  values  of  g}  for  negative  j  resulted 
in  a  65  element  kernel  for  use  in  Equation  15.  Figure  3  indicates  the 
shape  of  this  kernel.  The  horizontal  scale  is  noted  in  terms  of  time 
(tj  *  jAt).  The  simulated  linear  subsystem  appears  to  have  a  memory 
of  about  3  seconds.  Convolution  of  this  kernel  with  cosinusoidal  input 
and  harmonic  analysis  of  results  yielded  correspondence  with  the  analy¬ 
tical  simulation,  Equation  10,  to  about  3  significant  figures. 

Truncation  of  the  computed  values  of  the  quadratic  kernel,  gi,  , 

J  x 

was  also  carried  out  for  negative  values  of  the  indices.  This  resulted 
in  a  65  x  65  element  kernel  (j,k  »  0. ..64).  Convolution  of  this  kernel 
with  cosinusoidal  input  and  harmonic  anlaysis  of  results  also  yielded 
correspondence  to  three  significant  figures  with  the  analytical  simu¬ 
lation,  Equation  11,  for  G2(gj,w),  and  G^U.-w).  Figure  4  is  an  "iso¬ 
metric"  picture  of  the  significant  part  of  the  quadratic  kernel.  The 
intersections  of  the  lines  denote  the  actual  values  of  the  kernel, 
and  the  values  for  j  and  k  greater  than  49  have  been  omitted  from  the 
picture  because  they  are  very  small.  The  "horizontal"  scales  are 
given  in  terms  of  time  (tj  *  jAt  and  t  *  kAt),  the  point  0,0  being 
at  the  "near"  corner.  Though  the  horizontal  rotation  of  the  surface 
was  made  slightly  beyond  45  degrees  to  improve  the  appearance,  the 
right  to  left  symmetry  (exact  in  the  numerical  work)  is  apparent. 

The  quadratic  sub-system  also  has  a  memory  of  about  3  seconds. 


Because  the  third  term  of  the  discrete  time  domain  prediction 
Equation  17  involves  on  the  order  of  3N3  multiplications  per  response 
point  it  was  of  considerable  practical  importance  to  truncate  gj^ 
as  much  as  reasonably  possible.  Omitting  the  small  values  computed  for 


negative  values  of  j ,  k  and  4  results  in  a  65  x  65  x  65  matrix  for 
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Examination  of  the  computed  values  of  this  discrete  kernel  for  large 
values  of  the  indices  suggested  that  it  would  be  reasonable  to  truncate 
the  kernel  at  j,  k  and  l  of  49  so  that  the  final  truncated  kernel  became 
a  50  x  50  x  50  matrix  (j ,  k,  l  =  0...49).  Convolution  of  this  truncated 
kernel  with  cos i nu i soi da  1  input,  and  analysis  per  Equation  7b,  yielded 
estimates  of  G^  (oj  ,u> ,  -id  )  and  G3(ua,cu,w)  which  compared  with  the  analytical 
values  from  Equation  12  to  between  1  and  2  significant  figures.  This 
was  less  precision  than  was  desired  and  it  was  suspected  that  the  error 
was  due  to  the  truncation.  However,  a  65  x  65  x  65  kernel  requires  more 
than  twice  the  number  of  multiplications  in  Equation  17  as  a  50  x  50  x  50 
The  less  than  desired  precision  was  accepted  for  the  present  exploratory 
study  in  return  for  a  significant  economy. 

It  was  of  considerable  interest  to  see  what  a  cubic  impulse 
response  function  looks  like.  The  function,  gj^,  is  four-dimensional 
so  that  the  best  approach  appeared  to  be  to  look  at  the  three  dimensional 
functions  which  can  be  formed  if  one  of  the  three  time  delay  arguments 
is  held  constant.  Figure  5  shows  nine  of  the  50  possible  surfaces  which 
may  be  examined  by  holding  Z(t3)  in  gj^  constant.  The  quas  i  -  i  sometr  i  c 
projection  is  the  same  for  each  surface  and  is  almost  the  same  as  that 
used  in  Figure  4  for  the  picture  of  the  quadratic  frequency  response 
function.  The  "vertical"  scales  in  each  surface  are  the  same.  The 
"horizontal"  scales  in  each  case  are  annotated  in  terms  of  time  delays 
t1  and  t2 ,  and  the  constant  value  of  t3  is  noted  to  the  right  of  each 
surface.  The  surface  for  t3  =  0  is  essentially  flat  to  the  resolution 

shown,  as  is  that  for  t  =  3.0625  U  *  49).  Between  t,  =  2  and  3  the 

3  J 

"ripples"  shown  in  the  surface  at  the  top  of  the  figure  gradually  die 
out,  so  that  the  surfaces  shown  include  the  significant  parts  of  the 
cubic  impulse  response.  Since  by  assumption: 

g3(t2,t1,t3)  -  g3(t1,t2,t3) 

there  is  exact  right  to  left  symmetry  in  each  of  the  surfaces.  The  numer 

leal  results  show  the  required  symmetry  in  the  t3  direction  as  well,  so 

that  the  same  pictures  would  have  been  obtained  had  t  or  t  been  held 

1  2 

constant  and  the  surfaces  drawn  as  functions  of  the  remaining  two  time 
delays.  It  appears  that  the  cubic  sub-system  has  about  a  3  second  memory 
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Qualitatively,  there  appears  to  be  no  feature  of  any  of  the  surfaces 
shown  in  Figure  5  which  are  different  than  might  be  expected  in  a 
quadratic  impulse  response. 
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TIME  DOMAIN  COMPUTATIONS:  PRACTICAL  MATTERS 

As  is  easily  appreciated,  direct  implementation  of  the  time  domain 
computation  of  Equation  14  with  a  (50  x  50  x  50)  cubic  kernel  is  relatively 
expensive,  and  is  in  fact  wasteful  of  time  and  computer  memory  because  of 
the  symmetries  of  the  impulse  functions.  (The  convolutions  mentioned  in 
the  last  section  were  carried  out  with  the  methods  to  be  described.) 

There  is  nothing  which  can  be  done  to  speed  up  the  linear  part  of 

the  computation  (Y1 (n) ,  Equation  15)  and  in  practice  it  is  implemented 

as-is.  Some  re-organization,  and  the  two-fold  symmetry  of  the  quadratic 

frequency  response  function,  g-  allows  the  number  of  multiply-add 

J  k 

operations  in  Equation  16  to  be  reduced  by  a  factor  of  something  more 
than  two.  If  the  limits  of  the  j  and  k  summations  in  Equation  16  are 
defined  as: 

p  ^  j  <  r 
p  <  k  **  r 


the  Equation  16  may  be  re-written  as: 
Y2(n)  =  gzrr  X2(n  -  r) 


r-1 

+  l  X (n  -  j) 
j=P 


9h .  X(n  -  j )  +  V  2gt  X  (n 
JJ  k=j+1  JK 


-  k) 


(24) 


In  practice  it  is  advantageous  to  " tr iangu lar i ze"  the  quadratic  kernel 

for  use  in  Equation  24.  This  just  takes  the  form  of  creating  a  one 

dimensiona’  array  of  values  of  gT,  in  the  order  required  by  the  details 

J  k 

of  the  formula.  Equation  24,  and  absorbing  the  factor  of  two  in  each  of 
the  elements  wnere  j  /  k. 

It  is  clearly  most  important  to  take  advantage  of  the  six-fold 
symmetry  of  gj^  in  Equation  17: 

q  ^ 

y  j  ki 


9jik 


3 

kji 

3 

ki  j 


g  3 

9ijk 


=  g 


ikj 
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As  in  the  quadratic  case  the  limits  on  the  summations  in  Equa¬ 
tion  17  may  be  defined  as: 

p  <  j  ^  r 

p  <  k  <  r 

p  <  l  <  r 

With  these  conventions  Equation  17  may  be  re-written  as: 


Y3(n)  =  X2 (n-p) 

+  X (n-p- 1 ) 


gppp  X(n-P}  +  3gpp>p+1  X(n-p-l} 


x{n-p-l)  [g3+1>p+1)p+1  X(n-p-l)  +  3g3>p+1>p+1  X(n-p)_ 


+  l  X(n-j) 
j=p+2 


X(n-j)  jgjjj  X(n-j)  +  £  3g?jk  X(n-k)J 


+  J  3gjkk  X2(n_k) 

k=p 

j-1  k-1 

+  l  X(n-k)  l  6g?  X(n-£) 

k=p+1  JK* 


e=p 


(25) 


"Tr iangulari zation"  of  the  kernel  in  this  case  reduces  the  size  of  the 
kernel  by  a  factor  of  6  and  allows  the  factors  of  3  and  6  to  be  absorbed. 
The  scheme  is  the  same  as  for  the  quadratic  kernel,  a  one  dimensional 
array  of  values  of  gjk^  is  created  in  the  order  required  by  the  details 
of  the  formula,  Equation  25.  This  approach.  Equation  25,  reduces  the 
number  of  multiply-add  operations  per  response  point  required  for  Y3(n) 
with  a  50  x  50  x  50  kernel  from  375000  for  Equation  17  to  24600,  a 
reduction  by  a  factor  of  15. 


In  actual  programming  account  has  to  be  taken  of  the  memory  of 
the  system,  and  the  result  is  that  the  first  valid  point  of  the  response 
to  X(n),  n=1,2...  is  Y(l+r),  and  the  excitation  must  be  available  for 
p  points  past  the  last  response  point  computed.  In  the  present  instance 
only  the  valid  range  of  response  and  excitation  was  retained. 
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SIMULATION  OF  RESPONSE  TO  RANDOM  EXCITATION 

Within  the  objectives  of  the  present  work  it  was  as  important  to 
obtain  a  qualitative  idea  of  the  nature  of  cubic  response  to  random 
excitation  as  it  was  to  generate  samples  with  which  to  work  later  on. 
Accordingly,  given  the  kernels  of  the  simulated  system,  and  the  algorithms 
just  described,  it  remained  to  generate  random  excitation  sequences  and 
turn  on  the  computer. 

For  this  purpose  it  was  convenient  to  use  a  subroutine  left  over 

.L 

from  some  previous  work.  Reference  23  .  What  the  subroutine  does  is  to 
generate  digital  sequences  from  a  pseudo  random  band  passed  Gaussian  zero 
mean  process.  In  the  routine  utilized  the  numerical  band  pass  filter  is 
configured  so  that  the  resulting  sequence  nominally  has  the  spectral 
form: 

U>)  -  5cj2  ^  Exp{-1 ,25(u  /u)4}/oj5  (26) 

where  =  the  single  sided  spectrum 

=  spectrum  area 
=  variance 

wq  =  modal  frequency 

and  the  values  of  variance  and  modal  frequency  may  be  specified.  This 
spectral  form  is  the  same  as  those  of  the  ITTC  two  parameter,  Pierson- 
Moskowitz,  and  Bretschneider  wave  spectra.  The  algorithm  is  not  unusual. 
Random  numbers  from  a  uniform  distribution  are  first  generated  by  a 
standard  computer  system  utility.  These  numbers  are  conceptually  assumed 
to  represent  a  Gaussian  probability  and  a  numerical  approximation  is  then 
used  to  generate  zero  mean  Gaussian  deviates  with  unity  variance  and  a 
"white"  spectrum.  The  approach  has  been  checked  statistically  with  very 
large  samples  and  found  to  reasonably  represent  the  required  zero  mean 
white  Gaussian  process  within  5  standard  deviations.  The  next  step  is 
to  filter  the  sequence  numerically  to  shape  the  signal,  and  finally  to 
carry  out  time  wise  interpolation  to  generate  a  sequence  with  the  required 
sampling  interval. 


23.  Dalzell,  J.F.,  "A  Note  on  the  Distribution  of  Maxima  of  Ship  Rolling," 
Journal  of  Ship  Research,  Vol .  17,  No.  A,  December  1973- 
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In  the  present  case  the  sampling  interval  was  fixed  at  0.0625 
seconds  by  the  previous  derivation  of  the  impulse  responses.  It  seemed 
reasonable  for  a  first  attempt  to  set  =  2^  so  that  the  peak  of  the 
spectrum  and  the  linear  resonance  would  be  aligned.  For  present  purposes 
it  was  desirable  to  compute  response  for  various  excitation  levels,  and 
it  made  sense  to  take  account  of  the  homogeniety  property.  In  particular, 
if  Y1 (n) ,  Y2(n)  and  Y3(n)  {Equations  15  through  17)  are  computed  for 
some  nominal  level  of  excitation  and  stored,  then  the  response  to  an 
excitation  which  is  a  factor  "f"  times  the  nominal  is: 

fYl(n)  +  f2Y2(n)  +  f3Y3(n) 

which  is  to  say  that  once  the  computation  of  Equations  15  through  17  is 
made  for  one  excitation  level.  Equation  14,  the  total  response,  may  be 
evaluated  for  any  number  of  linearly  proportional  excitations  with 
trivial  extra  expense.  On  this  basis  it  was  decided  the  nominal  vari¬ 
ance,  cr2,  of  Equation  25  should  be  0.0625  so  that  o  is  0.25  and  the 

X  X 

"significant  height"  of  the  basic  excitation  would  be  unity. 

What  remained  to  be  settled  was  the  sample  duration.  The  routine 
used  had  provision  for  11  widely  separated  entry  points  into  the  under¬ 
lying  pseudo  random  number  sequence.  It  seemed  realistic  then  to  gener¬ 
ate  11  samples  of  a  duration  similar  to  that  ordinarily  achieved  in 
towing  tank  work.  Thus  2100  point  (131  second)  durations  were  chosen 
for  each  sample.  (This  choice  results  in  about  150  "wave"  encounters, 
and  broadly  corresponds  to  what  is  often  achieved  in  towing  tank  exper¬ 
iments  at  zero  ship  speed.) 

Figures  6a  and  b  indicate  the  time  histories  resulting  for 
Sample  1  with  the  basic  excitation  level  (o*  -  0.25).  The  entire 
sample  is  shown  in  two  figures  so  that  individual  flucuations  are  more 
visible.  At  the  top  of  each  is  the  excitation,  X(t).  The  three  frames 
in  the  middle  show  the  three  components  of  the  response,  Yl(t),  Y2 ( t ) 
and  Y3(t).  Finally,  at  the  bottom  the  sum,  or  total  response  Y(t), 
Equation  14,  is  shown.  It  should  be  noted  that  the  vertical  scales 
for  the  response  components  are  different;  the  scale  choice  was  made 
so  as  to  show  the  nature  of  each  component. 
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Qualitatively,  the  linear  and  quadratic  components  are  exactly 
as  expected.  The  linear  response  appears  statistically  symmetric  and 
contains  a  narrower  band  of  frequency  components  than  the  excitation. 

The  quadratic  response  is  quite  asymmetric,  includes  a  non  zero  mean, 
and  response  at  very  low  frequencies  as  well  as  a  visible  high  frequency 
response  which  is  predominately  twice  the  typical  linear  response 
frequency . 

No  ideas  about  the  qualitative  nature  of  the  cubic  component, 

Y3(t),  were  available,  and  it  was  this  question  which  occasioned  the 
plotting  of  time  histories.  As  is  clear  from  Figure  6,  the  simulated 
values  of  the  cubic  component  contain  "bursts"  of  response  which  appear 
symmetric  and  largely  to  have  the  same  frequency  content  as  the  linear 
response.  These  bursts  of  significant  cubic  response  appear  to  be 
associated  with  groups  of  waves  in  the  excitation,  and  with  significant 
excursions  of  the  quadratic  response,  Y2(t).  No  obvious  high  frequency 
components  appear  in  Y3(t).  In  retrospect,  some  of  this  qualitative 
behavior  of  the  cubic  component  of  response  might  possibly  have  been 
predicted.  The  simulated  cubic  frequency  response  function,  (ujj  ,00^  ,w3) , 
Equation  12,  contains  the  quadratic  frequency  response  function.  Thus 
some  qualtiative  relationship  should  result.  It  was  found  in  the  pre¬ 
vious  numerical  evaluations  that  the  cubic  response  at  the  third  harmonic 
of  excitation  frequency  was  an  order  of  magnitude  smaller  than  the  cubic 
response  at  the  excitation  frequency,  and  a  general  absence  of  high 
frequencies  from  the  cubic  response  could  have  been  expected  on  this 
bas i s . 

It  was  of  interest  to  see  if  the  qualitative  features  of  the  sim¬ 
ulated  nonlinear  response  held  from  sample  to  sample.  To  this  end  a 
number  of  the  remaining  samples  were  also  plotted.  Figures  7  and  8  are 
examples.  Sample  4  was  found  to  contain  the  highest  excursion  of 
cubic  response  of  all  the  eleven  samples.  This  largest  excursion  of 
Y3(t)  is  shown  at  about  t=5  seconds  in  Figure  7.  Sample  11  was  found 
generally  to  have  the  lowest  level  of  cubic  response,  and  half  of  this 
sample  is  shown  in  Figure  8.  The  qualitative  nature  of  the  simulated 
cubic  response  to  random  excitation  appears  to  hold  throughtout  the 
simulated  data. 
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A  general  comparison  of  the  total  response,  Y(t),  in  Figures  6 
through  8  with  the  corresponding  linear  component,  Y 1  ( t ) ,  indicates 
almost  no  qualitative  difference.  At  the  basic  excitation  level  the 
nonlinearities  are  too  small  to  have  much  effect.  Those  quantitative 
changes  which  are  visible  are  most  easily  seen  in  the  first  part  of 
Figure  7.  There  the  largest  nonlinear  responses  in  the  entire  simu¬ 
lation  produce  a  non  symmetric  response  and  perhaps  slightly  increase 
the  amplitude  of  the  fluctuations  relative  to  those  of  the  linear 
component . 

It  was  also  of  interest  to  see  the  qualitative  effect  "oon  the 
total  response,  Y ( t ) ,  of  an  increase  in  the  excitation  leve.,  c  . 

Figures  9a, b,  10  and  11  indicate  this  for  an  excitation  four  times  the 
nominal  (a^  =  1.0)  for  the  cases  shown  in  Figures  6  through  8.  In 
accordance  with  the  previous  discussion,  the  synthesis  of  these  results 
amounts  to  multiplication  of  the  basic  excitation  and  linear  response 
by  4,  multiplication  of  quadratic  response  by  16,  and  multiplication 
of  cubic  response  by  64  before  making  the  summation  for  Y(t).  Thus 
the  time  histories  in  the  top  four  frames  of  Figure  9a  for  instance 
are  exactly  the  same  as  those  in  Figure  6a  except  for  a  scale  change. 

The  qualitative  difference  is  in  the  bottom  frame,  Y(t).  It  is  clear 
from  Figures  9  through  11  that  once  the  excitation  level  is  high  enough 
there  are  significant  qualitative  differences  between  the  linear  and  the 
the  nonlinear  response  to  randan  excitation,  and  the  differences  center 
about  occasional  non  symmetric  large  amplitude  excursions  which  change 
the  general  appearance  of  the  time  history. 

Though  some  of  the  features  of  the  total  simulated  response  in 
Figures  9  through  11  may  be  a  bit  more  extreme  than  some  response  his¬ 
tories  obtained  in  towing  tank  experiments,  some  of  the  same  behavior 
is  occasionally  seen. 
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THE  SPECTRUM  OF  RESPONSE  TO  RANDOM  EXCITATION 


Regardless  of  the  possibility  of  non  1 i near i t i es ,  the  scalar  spectra 
of  excitation  and  response  are  the  fundamental  state  of  art  approach  to 
the  characteri zation  of  the  random  case.  The  theoretical  form  of  the 
spectrum  for  the  present  model.  Equation  3,  may  be  written  down  from 
References  A  and  7.  For  subsequent  purposes  it  is  convenient  tc  cite 
the  results  firc,t  in  a  particular  two-sided  form.  If, 


S  (w)  =  excitation  spectrum 

XX 

S  (w)  =  response  spectrum 

yy 

where  both  are  double  sided  and  normalized  so  that  their  integral  over 
frequency  is  equal  to  (2m)  times  variance,  then  the  response  spectrum 
for  the  1 inear+quadra t i c+cub ic  mode!  may  be  written  in  terms  of  frequency 
response  functions  as  follows: 

Syy(u)  =  S2U)  +  S,(w)  +  S3(u)  (27) 

and : 

S.(w)  =  SvU)  |G  .(ui)  +  y-  G  (w,v,-v)  S  (v)  dv  |  -  (28) 

*  XX  ^  j  j  XX 


s2  (ui) 


TT 


(lD~U,u) 


I  O 

I 


S  (uj-u) 
xx 


S  (u)du 

XX 


(29) 


S3(w) 


G  (uj-v-u , v ,u)  | : 


S  (w-v-u) 
xx 


S  (v)  S  (u)dudv 
xx  xx 


(30) 


8y  way  of  comment,  if  the  cubic  frequency  response  function  is  zero 
the  form  of  the  spectrum  is  the  same  as  that  for  the  linear  plus  quadratic 
system  of  Reference  15.  If  both  cubic  and  quadratic  frequency  response 
functions  are  zero,  the  result  is  the  conventional  linear  estimate. 

In  practice,  as  opposed  to  theoretical  manipulations,  the  one 
sided  realizable  spectrum  with  area  equal  to  variance  is  required,  and 
the  above  equations  may  be  written  in  terms  of  one  sided  spectra  by 
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rep  lac i ng 


S  (■*,)  and  S  ) 
xx  yy 

U  (uj)  =  the 

xx 

U  (w)  =  the 

yy 


by  -U  U)  and  ~U  (-)  where 
xx  yy 

single  sided  excitation  spectrum 
single  sided  response  spectrum 


Carrying  this  and  some  simple  variable  changes  out  results  in  expressions 
for  the  single  sided  spectra  which  are  analogous  to  Equations  27  through  30: 


Uyy(u)  =  U1  (w)  +  U2  (oj)  +  U3(ui) 


(3D 


where  (w  >  0)  and: 


U1  (w)  =  (uj  )  (G^w)  +  3 


G,(w,v,-v)  U  (v)dv| 2 


xx 


(32) 


U2(w)  = 


G^ {  (w-v) /2,  (oj+v) /2}  ]  2  l>^(  |  uj-v  | /2)  ( | u+v  | /2) dv  (33) 


U3(«o)  =  1.5 


JJ 


|G,(oj-v-u,v,u)  l2  Uvv(  |  ou-v-u  1 )  Uvv(jv|)  Uvv(|u|)dudv 

(3M 


Implicit  in  Equation  3^  is  a  failure  to  find  a  change  of  variables 
which  would  allow  reduction  of  the  range  of  integration.  Given  the 
excitation  spectrum  and  expressions  for  the  frequency  response  functions. 
Equations  32  through  3^  are  computable,  though  Equation  3^  is  tedious. 

If,  as  in  the  present  case,  there  is  interest  in  response  spectra  for 
excitation  which  is  linearly  proportional  to  some  nominal  excitation, 
the  integrations  in  Equations  32,  33  and  3^  need  only  be  done  once. 

(Once  the  integrations  are  performed  for  a  given  U^x(w),  the  integrals 
for  an  excitation  a  factor  "f"  times  the  original  are  the  initial  results 
times  (f2lc)  where  k  is  the  number  of  times  llxx(u))  appears  in  the  integrand.) 

It  was  of  interest  both  to  see  if  all  the  above  works  and  to 
obtain  some  idea  of  the  relative  influence  of  the  various  terms.  Thus 
it  was  decided  to  compose  simulated  time  histories  for  various  levels 
of  excitation,  carry  out  conventional  spectral  analyses,  and  finally 
do  the  computations  implied  in  Equations  31  through  31*  for  comparison. 
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The  first  step  in  this  operation  was  to  produce  time  histories  of 
X(t)  and  Y(t)  for  excitations  1/2,  1,  2,  3  and  4  times  the  nominal 
excitation  specified  in  the  last  section.  As  explained  previously, 
given  the  stored  time  history  components  Yl(t),  Y2(t)  and  Y 3 ( t )  for 
the  nominal  excitation,  the  production  of  X(t)  and  Y(t)  for  each  excit¬ 
ation  level  is  a  trivia)  exercise.  Once  it  is  done,  the  resulting  time 
series  appear  more  or  less  like  any  samples  of  random  input  and  output, 
and  conventional  data  reduction  techniques  can  be  employed. 

Spectrum  analysis  of  each  time  history  was  first  carried  out  by 
the  conventional  Fast  Fourier  Transform  frequency  smoothing  technique 
noted  in  Reference  24  .  Each  time  series  was  truncated  from  2100  to 
2048  points,  corrected  to  zero  mean,  tapered  with  the  Tukey  10%  cosine 
taper  and  directly  transformed  with  the  FFT.  The  raw  spectral  estimates 
resulting  were  averaged  in  groups  of  13  at  intervals  of  7  to  result  in 
spectral  estimates  with  26  degrees  of  freedom  and  a  frequency  resolution 
of  0.344  rad/sec.  For  each  excitation  level  there  were  11  such  spectra 
resulting  from  this  operation  on  excitation  and  response.  Since  each 
time  history  is  a  sample  from  a  stationary  random  process  it  was  con¬ 
sidered  reasonable  to  also  carry  out  ensemble  smoothing  over  the  11 
samples  for  each  excitation.  Thus  the  final  result  for  each  of  the 
five  excitations  was  a  single  smoothed  "observed"  excitation  spectrum 
and  a  single  "observed"  response  spectrum.  Each  of  these  spectra  had 
286  degrees  of  freedom  per  spectral  estimate,  which  implies  90%  con¬ 
fidence  bounds  on  the  estimates  of  +15%  and  -12%.  Total  degrees  of 
freedom  for  the  excitation  spectra  were  in  excess  of  3000,  which  implies 
90%  confidence  bounds  on  variance  of  ±4%,  and  in  fact  the  estimated 
variance  was  within  2%  of  the  theoretical  value  specified  in  the  original 
time  domain  simulation.  Total  degrees  of  freedom  for  the  "observed" 
response  varied  between  1500  and  2500,  which  would  result  in  confidence 
bounds  on  variance  of  +5%  if  the  response  could  be  considered  Gaussian. 
Table  4  summarizes  the  "observed"  standard  deviations  from  the  data  re¬ 
duction  procedure.  It  should  be  noted  that  the  "observed"  values  of 


A 

24.  Bendat,  J.S.  and  Piersol,  A.G.,  "Random  Data:  Analysis  and  Measure¬ 
ment  Procedures,"  John  Wiley  6  Sons,  1971. 
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TABLE  4 

"OBSERVED"  STANDARD  DEVIATIONS  OF  EXCITATION 
AND  RESPONSE  FOR  VARIOUS  NOMINAL  EXCITATION  LEVELS 


Factor  on 

Nominal  Excitation 

0.5 

1 .0 

2.0 

3.0 

4.0 

Nom  i  na  1  a 

X 

0.125 

0.25 

0.50 

0.75 

1 .00 

"Observed"  a 

X 

0.124 

0.248 

0.496 

0.744 

0.992 

"Observed"  a 

y 

0.181 

0.356 

0.673 

0.971 

1  .396 

a  /a 

1.46 

1.44 

1.36 

1.30 

1  .41 

y  * 
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standard  deviation,  o  ,  march  upward  exactly  as  expected--the  operations 
carried  out  on  the  excitation  time  histories  can  only  result  in  an  i  n- 
crease  in  standard  deviation  proportional  to  the  assumed  factor  on 
nominal  excitation.  The  "observed"  response  standard  deviations  are 
not  exactly  proportional  to  the  assumed  factor  on  excitation  as  must 
be  expected  for  nonlinear  systems. 

The  next  step  was  to  carry  out  the  "prediction"  operations  of 
Equation  31  through  3^.  The  frequency  response  functions  for  the  sim¬ 
ulated  system.  Equations  10  through  12,  had  been  programmed  for  earlier 
operations  to  that  numerical  values  for  the  functions  could  be  obtained 
easily.  Apart  from  carrying  out  the  operat ions,  the  only  other  question 
was  whether  to  use  the  theoretical  expression  for  U^Cu),  Equation  26, 
or  the  "observed"  excitation  spectrum  from  the  data  reduction  procedure. 

It  was  guessed  that  use  of  the  "observed"  excitation  spectrum  might 
wash-out  some  of  the  residual  effects  of  statistical  variability  and 
the  results  of  the  data  reduction  were  used  as  the  excitation  in  the 
evaluation  of  Equations  31  through  3^  for  the  five  cases  of  interest. 

Figures  12  through  16  summarize  the  results  of  the  comparison. 

Each  figure  has  four  frames  with  frequency  the  abcissa  and  spectral 
density  the  ordinate.  The  "observed"  excitation  spectrum  is  shown  in 
the  uppermost  frame,  and  the  "observed"  response  spectrum  in  the  lower 
most,  both  plotted  as  boxed  points  connected  by  straight  lines.  The 
"predicted"  response  spectrum,  Equation  31,  is  overplotted  in  the 
lower  frame.  The  middle  frames  indicate  the  U1 (u) ,  U2(w)  aid  U3(w) 
components  of  the  prediction.  Equations  32  through  3^.  Finally,  as 
an  aid  in  interpretation,  the  purely  linear  estimate  of  the  spectrum 
is  indicated.  This  is  just  what  would  be  obtained  if  both  quadratic 
and  cubic  nonlinearities  were  ignored,  or: 

Uxx^)  1 G  j  (w)  | ' 

In  Figure  12,  for  the  lowest  excitation  level,  it  is  clear  that 
the  nonlinearities  would  be  of  no  practical  importance.  U2 (-o )  and  U3(w) 
are  invisible  and  the  influence  of  the  cubic  nonlinearity  upon  111 (w)  is 
very  slight.  The  agreement  between  observed  and  predicted  response  spectra 
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FIGURE  14  OBSERVED  AND  PREDICTED  SPECTRA  OF  RESPONSE 


OF  THE  SIMULATED  SYSTEM:  NOMINAL  o 
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is  almost  embarrassingly  good.  Very  much  the  same  situation  appears  in 
Figure  13  for  the  nominal  excitation,  r ^  =  0.25.  The  very  low  frequency 
response  spectral  density  due  to  quadratic  i nteract ions ,  U2  D  ,  just 
starts  to  appear,  and  the  influence  of  the  cubic  nonlinearity  on  U1 (~) 
is  slightly  greater  than  in  Figure  12. 

In  Figure  1**  the  influence  of  the  U3(w)  component,  Equation  3^, 
just  starts  to  appear.  Figures  15  and  16  for  =  0.75  and  1.0  indicate 
the  form  of  the  contributions  to  the  spectrum  of  all  three  components. 
The  nonlinear  part  of  U1 (o) ,  Equation  32,  can  significantly  change  what 
might  be  expected  from  purely  linear  considerations.  U2(a,),  the  quad¬ 
ratic  part,  make-  a  contribution  to  very  low  frequency  response, 
perhaps  to  frequencies  in  the  excitation  range,  and  contributes  some 
high  frequency  response  at  frequencies  about  a  factor  of  2  above  the 
excitation  peak.  U3(w),  the  pure  cubic  part,  is  only  significant  at 
the  highest  excitation  level,  and  makes  a  contribution  in  the  range 
of  the  peak  of  the  spectrum. 


It  is  impossible  to  generalize,  but  if  the  simulated  system  is 
representative,  it  may  be  that  in  many  practical  seakeeping  problems 
only  the  first  component  of  the  spectrum,  U1 (w) ,  will  be  of  importance. 


One  of  the  most  common  operations  in  analysis  of  random  wave 
experiments  is  to  estimate  the  modulus  of  the  linear  frequency  response 
function  by  taking  the  square  root  of  the  ratio  of  observed  response 
and  excitation  spectra: 


|Ga  (o>) 


estimated 


U  (w)/U  U)  3 
yy  xx 


This  operation  was  performed  using  the  "observed"  response  and  excitation 
spectra  just  presented.  The  results  are  shown  in  Figure  17-  In  accor¬ 
dance  with  usual  practice  no  results  are  given  for  frequencies  where 
the  excitation  spectrum  was  less  than  10%  of  its  peak. 

The  equivalent  linear  response  appears  to  vary  systematically 
with  excitation  level  o^.  As  excitation  level  increases  the  results 
deviate  more  and  more  from  the  form  of  the  simulated  linear  response 
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function,  until  finally  there  is  no  resemblence.  The  changes  at  the 
ends  of  the  frequency  range  shown  in  Figure  17  are  clearly  due  to  quad¬ 
ratic  response  (U2(^)).  The  changes  in  midrange  resemble  those  noted 
in  Figure  1  for  the  influence  of  the  cubic  nonlinearity  on  the  response 
to  single  tone  excitation. 

As  in  the  discussion  of  the  qualitative  features  of  the  response 

to  random  excitation  of  the  simulated  cubic  system,  the  response  spectra 

for  excitations  of  o  =  0.75  and  1.0  deviate  from  linear  expectations  by 

x 

amounts  which  might  be  considered  extreme  by  comparison  with  towing 
tank  experiments.  However,  the  magnitudes  of  nonlinear  distortion  for 
values  of  up  to  0.5  is  in  line  with  what  has  been  occasionally  ex¬ 
perienced.  It  thus  appeared  that  the  simulation  was  reasonably  realistic. 
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IDENTIFICATION  THEORY  FOR  THE  RANDOM  CASE 


It  is  clear  from  the  preceding  section  that  the  ordinary  scalar 
spectrum  analyses  of  random  excitation  and  response  provides  no  means 
of  i dent i f i cat  ion  of  frequency  response  functions.  What  scalar  spectrum 
analysis  can  provide  is  an  indication  if  non  1 i near i t i es  are  serious  or 
not  in  the  case  where  data  is  available  for  a  range  of  excitation  levels. 
The  problem  of  interest  here  is  whether  or  not  the  three  postulated 
frequency  response  functions  arising  from  the  basic  model,  Equation  5, 
may  be  identified  from  samples  of  random  excitation  and  response. 

To  begin,  the  assumptions  about  the  excitation  previously  noted 
should  be  re-stated.  It  is  accepted  that  the  excitation  is  strictly 
stationary  with  bounded  moments  of  all  orders.  For  present  purposes 
the  conventional  further  assumptions  about  a  wave  process  will  also 
be  made;  that  is,  the  excitation  will  be  assumed  to  be  an  ergodic 
Gaussian  zero  mean  process.  The  autocorrelation  function  of  the 
excitation  will  be  denoted: 


Rxx(t)  «  X(  t)  X(t  -  T) 


(35) 


where  the  overbar  denotes  the  temporal  mean,  and  by  assumption  the 
statistical  expectation.  The  two  sided  spectrum  of  the  excitation  will 
be  defined  as: 


S 

xx 


(<*>) 


R  (t)  Exp  (-iwx)dT 


(36) 


where,  as  before,  the  integral  of  the  spectrum  is  2 tt  times  the  variance, 
Rxx(0).  Given  the  Gaussian  zero  mean  assumption,  all  the  expectations 
of  higher  order  products  are  either  zero  or  are  defined  in  terms  of 
Rxx(’r)»  Equation  35.  Appendix  A  summarizes  the  expectations  of  up  to 
sixth  order  products  from  Reference  25  . 


The  general  analytical  approach  in  the  present  section  of  the 
report  will  be  to  hypothesize  certain  statistical  time  domain  moments 


7C 

25.  Laning,  J.H.,  and  Battin,  R.H.,  "Random  Processes  in  Automatic 
Control",  McGraw-Hill,  1956. 
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of  response  and  excitation,  and  then  invert  the  resulting  expressions 
to  the  frequency  domain.  In  this  operation  frequent  use  is  made  of 
an  n-d imens i ona 1  form  of  Parseval's  formula  as  given  in  Reference  k. 

For  convenience  this  formula  is  given  in  Appendix  B. 

The  first  response  moment  of  interest  of  the  cubic  model.  Equa¬ 
tion  5,  is  just  the  temporal  mean.  Taking  the  temporal  mean  of 
Equation  S' 

—r-  ' 

V ( t)  =  '  g.(t,)  X ( t  -  t  )dt, 

J  -  A  I  * 

r  f 

+  jj  9 2 ( t j » t_ )  Rxx(t,  -  t.)dt,dt. 

1 1  Rxx( -  t2)dt1dta 

=  2~  [  G7(uj,-ui)  Sxx(a;)du)  (37) 

where  the  derivation  follows  that  of  References  8,  14,  and  15.  It  may  be 
noted  that  this  result  for  the  cubic  model  is  exactly  the  same  as  that 
for  the  quadratic  model.  Reference  15,  because  the  odd  order  temporal 
means  are  zero  for  the  assumed  Gaussian  excitation  process. 

Because  cross-spectral  analysis  is  a  very  common  technique  in 
the  analysis  of  linear  random  processes,  it  is  of  interest  to  look  at 
the  cross  spectrum  between  excitation  and  the  nonlinear  response  model. 
The  two  sided  cross  spectrum  may  be  defined  as: 

S  (w)  =  [  R  (t)  Exp  (-iwT)dT  (38) 

yx  j  yx 

where  the  cross  correlation  function  between  X ( t )  and  Y(t)  is  defined: 

R  (t)  =  f(Y(t)  -  YTtT)  X(t  +  x)l  l39) 

Yx  L_  _i 

Substituting  the  model  for  Y(t),  Equation  5,  into  Equation  39,  elimin¬ 
ating  temporal  means  of  odd  products  of  Gaussian  variables,  applv:ng 


1(t1,t2,t3)  x"Ct  -  tj)  X  ( t  -  tj  X  ( t  -  t,’5dt1dt;dt3 
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the  expression  for  fourth  order  expected  values,  Appendix  A,  and  taking 
advantage  of  the  symmetry  of  the  impulse  functions  results  in  a  further 
expression  for  the  cross  correlation: 


VT)  '  ..  Rxx<:  * 


♦  3  j:j  Rxx(t, 


tj  R  (t  +  tjdt.dt.dt, 

XX  3  1  -  3 

(AO) 


Now  substituting  Equation  AO  into  Equation  38,  and  applying  the  n- 
dimensional  Parseval  formula,  Appendix  B,  results  in  the  following 
ression  for  the  cross  spectrum: 


S 


yx 


03 


(u.) 


G3(w,v,-v)  Sxx(v)dv 


exp- 


(Al) 


Just  as  in  the  1 inear+quadrati c  case,  Reference  lA,  the  quadratic  con¬ 
tributions  drop  out.  The  result  in  Equation  Al  is  similar  to  the  first 
(and  possibly  most  important)  component  of  the  scalar  spectrum,  S^u), 
Equation  28. 

The  ordinary  coherency  between  excitation  and  response  is  defined 
as: 


Y2  (w)  =  |S  (uj)|2/{S  (u>)  S  (u>)}  (A2) 

yx  1  yx  1  yy  xx 

Now  substituting  Equation  27  and  Equation  Al  into  Equation  A2,  and 
noting  Equation  28,  the  theoretical  expression  for  the  coherence  becomes 


Y  2  .( u> ) 

yx 


s,(w)  +  nz 7 
1  +  ~  s^rrr — 


(A3) 


where  S^w),  S0(oj)  ind  S?(^)  are  defined  in  Equation  28  through  30. 
Clearly,  the  theoretical  coherence  cannot  be  exactly  unity  as  in  the 
fully  linear  case.  However  it  will  approach  unity  as  the  excitation 
level  decreases.  In  the  discussion  of  the  spectra  of  the  simulated 
process  of  the  last  section  it  was  noted  that  at  intermediate  levels  of 


» 
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excitation  S„(u.)  and  S,(u.)  may  be  negligible  in  comparison  with  (_) 
while  at  the  same  time  S^(^)  does  not  correspond  to  the  linear 
(negligible  excitation)  case.  Relative  to  the  interpretation  of 
actual  data,  this  suggests  that  there  may  be  situations  where  the 
sample  coherency  may  appear  reasonable  at  many  excitation  levels 
(suggesting  a  linear  response),  yet  the  apparent  frequency  response, 

S  (w)/Sxx(u)  ,  varies  with  excitation  level  (suggesting  a  nonlinear 
response).  Something  like  this  can  be  inferred  from  some  of  the 
results  in  Reference  26  ,  though  effects  of  short  samples  therein 
makes  the  conclusion  tenuous. 

It  is  clear  from  the  theoretical  form  of  the  cross  spectrum. 
Equation  41,  that  if  the  system  of  interest  has  a  significant  cubic 
response,  the  identification  of  the  linear  frequency  response  function, 
Gj(w),  by  conventional  cross  spectrum  analyses  would  require  a  very 
large  number  of  independent  experiments  at  various  excitation  levels. 
For  example,  if  it  was  thought  reasonable  to  discretize  the  integral 
in  Equation  41 ,  and  if  $xx(w)  could  be  represented  by  20  discrete 
values.  Equation  41  might  be  written  as  an  algebraic  equation  in  21 
complex  unknowns  for  each  value  of  w,  and  reasonable  results  might 
perhaps  be  obtained  with  data  from  40+  independent  experiments. 

Clearly,  Equation  41  is  unsatisfactory  as  a  basis  for  identification 
of  the  linear  function  unless  an  approach  for  G3(w,v,-v)  can  be  found. 

In  previous  work  with  the  1 inear+quadrat ic  model  an  identifi¬ 
cation  technique  called  cross-b i -spectra  1  analysis  was  used  to  identify 
quadratic  frequency  response  functions.  The  cross-b i -spectrum  may  be 
def i ned  as  fol lows : 

Ryxx(t4’t5)  EXP  ('iW!t4  '  ia,2t5)dt4dt5  (M° 

where  the  third  order  correlation  function  R  (t,  , t ,- )  is  written: 

yxx  4  •  5' 

Ryxx(t4>t5)  =  {Y(t)  -  YltT)  X(t  -  t4)  X(t  -  t5)  (45) 


J. 

26.  Dalzell,  J.F.,  "Some  Further  Experiments  on  the  Application  of  Linear 
Superposition  Techniques  to  the  Responses  of  a  Destroyer  Model  in 
Extreme  Irregular  Long-Crested  Head  Seas,"  Davidson  Laboratory 
Report  918,  September  1962. 


S  (w,  ,oi_ )  = 
yxx  1  2' 
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It  should  be  noted  that  this  definition  corresponds  to  that  of  Refer¬ 
ence  10,  and  that  the  definitions  given  in  References  14  and  15  are 
for  a  "modified"  cross-bi-spectrum  which  is  defined  in  terms  of  sum 
and  difference  frequencies  (  Wj  +  w.)  and  (wj  -  w  ) .  (The  "modified" 
cross-b i -spectrum  has  some  computational  advantage.) 

When  the  basic  model  for  the  response,  Y ( t ) ,  is  substituted  into 
Equation  45  there  results: 


g^tj)  x(t  -  t^)  x(t  -  t4)  xTt  -  t5)dt1 


g,(t1,t2)  X(t  -  tj)  X ( t  -  t2)  X ( t  -  t4)  X(t  -  t5)dt,dt2 


rr 


93(tl-t2-t3>- 


x(t  -  t\)  x (t  -  tj  x ( t  -  t3)  x ( t  -  t4)  x ( t  -  t'5T 


TTETR^It,  -  t5> 


dt  dt  dt 

12  3 


(46) 


Under  the  assumptions  given  for  the  excitation  the  expectations  of 
products  of  an  odd  number  of  excitation  variables  are  zero,  Appendix  A. 
Thus  the  terms  involving  linear  and  cubic  impulse  response  functions 
drop  out,  and  the  quadratic  impulse  response  is  isolated. 

Applying  the  expressions  for  the  fourth  order  expectation, 
Appendix  A,  and  noting  Equation  37,  Equation  46  becomes: 


R  (t,  ,tj  -  2 
yxx  4’  5 


g  (t.t)  R  (t  -  t,  )  R  (t  -  t Jdt ,dt 

2  1  2  XX  1  4  XX  2  5  12 


(47) 


Finally,  applying  the  Parseval  formula,  Appendix  B  to  Equation  47,  and 
noting  the  defintion  of  the  cross-b i -spectrum ,  Equation  44: 

S  (w  tw  )  -  26  (u  ,«  )  S  (u,  )  S  (w  )  (48) 

yxx  1’  2  2  l*  2  xx  l  xx  2 

From  this  result  an  estimator  for  the  quadratic  frequency  response 
function  becomes: 

Mv*2)  -  2lYl")TL) 

xx  1  xx  2 
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This  estimator  is  analogous  to  that  of  References  14  and  15-  The  impor 
tant  point  relative  to  the  present  investigation  is  that,  with  the 
Gaussian  assumption  on  the  excitation,  cross-b i -spectra  1  analysis  iso¬ 
lates  the  quadratic  frequency  response  whether  or  not  a  cubic  non- 
1 i near i ty  i s  present. 

Now  if  a  third  order  correlation  function  R  ?  (t,,t„)  is  formed 

xzxx  *  2 

by  replacing  the  response  in  Equation  45  by  the  square  of  the 
exci tation: 


R  2 

xzxx 


(tx,t2)  =  (x2(t)  -  x^TtT}  x ( t  -  t:)  x(t  -  t2) 


(50) 


and  the  double  Fourier  transform  is  taken  as  in  Equation  44,  there 
results  the  cross-bi-spectrum  between  excitation  and  excitation 
squared : 


S  o  (a)-.  )  = 

x2xx  1  2 


Rx2xx^tl,t2^  Exp (- itOj tx  -  iw2t2)dt1dt2 


-  2  5xx("l>  SxxK>  <S’> 

where  the  derivation  follows  that  in  Reference  15.  Accordingly  an 
alternate  identification  technique  may  be  written: 


;  ,  ,  _  ^xxS’1^ 

G2  (u)i  »w2)  =  s  2  C  ■.■o.TT 

X^XX  1*  2 


(52) 


Thus  the  quadratic  frequency  response  function  may  in  principle  be 
estimated  by  the  ratio  of  cross-b i -spectra ,  and  this  approach  was  noted 
in  Reference  15  as  resulting  in  improved  estimates. 


The  theoretical  approaches  to  the  reduction  of  random  data  which 
have  just  been  discussed  have  been  available  for  some  time.  The  next 
logical  step  in  the  present  work  was  to  look  for  an  approach  with  which 
there  might  be  hope  of  isolating  the  cubic  frequency  response  function, 
or  at  least  parts  of  it.  It  seemed  reasonable  to  proceed  by  analogy 
with  cross  and  cross-b i -spectra  1  theory  and  postulate  a  fourth  order 
correlation  function  of  the  form: 


Y ( t)  X(t  -  tu)  X(t  -  t5)  X(t  -  tj 
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It  should  be  clear  from  a  comparison  of  this  form  with  Equations  A5 
and  A6  that  the  quadratic  contributions  to  such  a  correlation  will  drop 
out  and  the  linear  and  cubic  parts  will  stay.  Substituting  the  model. 
Equation  5,  into  the  postulated  form,  expanding  the  expected  values  in 
accordance  with  Appendix  A,  and  taking  advantage  of  symmetry  results 
i n  the  fol lowi ng : 


Y(t)  X(t  -  tj  X ( t  -  t,)  X ( t  - 


I  ( t  )  [r  (t  -  t  )  R  (t  -  t  ) 

1  1  Lxx  1  4  XXV  5  6 

+  R  (t  -  t  )  R  (t,  -  t  ) 

XX  1  5  XX  4  6 

+  R  (t,  -  t,)  R  (t.  -  t J  dt, 

XX  1  s'  XX  4  5'J  1 


.11  VWS*  [Rxx(t, 


+  R  (t, 
xx  1 

+  R  (t„ 

xx  2 


+  R  (t.  - 
xx  1 

+  R  (t,  - 
xx  1 

+  R  (t  - 

xx  2 

+  R  (t.  - 
xx  1 

+  Rxx^l 


-  tZ)  Rxx(t3  - 

-  l3>  Rxxlt2  - 

-  l3>  Rxx(t3  * 


"xx<'3 
[3*  Rxx<£2 
C3>  Rxx(t3 
l2>  Rxx(t3 
Rxxlt2 


+  R  (t, 

xx  2 


-  t.)  R  (t,  - 
3'  xx  1 


t  )  R  (t 
4  xx  5 

t,  )  R  (tc 
4  XX  5 

t.  )  R  (tc 

4  XX  5 

‘s’  Rxx"x 

Rxx(t- 
tj  R  (t. 

5  XX  4 

tc)  R  (t, 

6  XX  4 

tj  R  (t. 

6  XX  4 

t,)  R  (t. 

6  XX  4 


t5>]dtldt2dt3 


+  6  g,(t, ,t0,t,)  R  (t-t,)R  (t-t)R  (t-t)dtdtdt 

JJJ  3  1*  2*  3  XX  1  4  XX  2  5  XX  3  6  1  2  3 

(53) 

The  problem  with  the  expression,  Equation  53,  is  that  it  contains  a 
linear  contribution  just  as  in  the  cross  spectrum,  Equation  40.  An 
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apparently  workable  approach  to  this  problem  is  provided  by  Reference  4 
in  a  treatment  of  partially  orthonormal  functional  expansions.  By  this 
means  a  fourth  order  correlation  function  was  postulated  as  follows: 

Ryxxx^t4,t5,t6*  =  '  'X(t  *  ‘  t  =  ^  X^  ‘ 


-  X  ( t  -  tj  X  ( t  -  tc)  X  ( t  -  tj 

4  z  5 


-  X(t  -  t5)  X  ( t  -  t4)  X  ( t  -  tj 

_ _  ~t  1 

-  x(t  -  t6)  X(t  -  t4J  X(t  -  t5T| j  (54) 

where  the  expection  operator  implies  that  temporal  means  are  to  be  taken 
Now  expanding  Equation  54  and  eliminating  terms  involving  the  expecta- 
t  ions: 


X(t  -  tTT  and  X{ t  -  t4)  X(t  -  t$)  X(t  - 

according  to  the  Gaussian  zero  mean  assumption: 


Ryxxx(t4>t5>t6)  =  X  ( t  -  t4)  X(t  -  t,.)  X(t  -  tj 


-  «xx(t4  -  t6)  v(t)  x(t  -  t5) 


-  -  te)  vTFrxrr^i 


Rxx(t.  -  t5)  Y(t)  X(t  -  t6: 


Mow  considering  terms  of  the  form  of  the  last  three  terms  of  Equation  55 


Rxx(tj  -  'k>  Vlt)  Xlt  -  h> 

'  |  W  Rxx(tl  -  \)  Rxx('j  ■  tk,dtl 

*  III  [Rxx<'l  -  Rxx<‘,  ■  <■>)  Rxx«j  -  ‘k> 

♦  "xx*1!  -  *3>  Rxx(t:  -  'i>  \x(tj  -  'k! 

*  Rxx(t:  •  '3J  Rxx(t.  -  ‘t>  “xx(tj  -  'k>]dt; 


dt.dt. 
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where  as  before,  the  model,  Equation  5,  has  been  substituted  and  the 
expected  values  expanded.  Comparing  Equations  55  and  56  with  Equa¬ 
tion  53  it  may  be  seen  that  the  effect  of  the  last  three  terms  of 
Equation  55  is  to  cancel  out  the  first  two  integrals  in  the  expression 


for  Y ( t)  X(t  -  t  )  X { t  -  t  )  X(t  -  t  ),  Equation  53.  Thus  the  postu- 

4  5  6 

lated  fourth  order  correlation  function  comes  down  to: 


R  ( t,  ,  t  ,  t  )  =6 
yxxx  4’  5  6 


93(tl*t-t3)' 


RXX(tl  -  t4)  Rxx(t2  “  t5)  Rxx(t3  '  t6)dtldt2dt3 

(57) 


Now  applying  the  Parseva!  formula,  Appendix  B,  to  Equation  57: 


R  (t,  ,tc,tj  ■ 
yxxx  4*  5’  6 


T27F 


G3(w1,“2*“3)  Sxx(w2)  SxxU3}  * 

Exp{i(u)1tu  +  w2t5  +  w3  tg)  }duldai7du3 


(58) 


Equation  58  is  in  the  form  of  a  triple  Fourier  transform  of  a  function 
of  frequency.  Inverting,  the  cross-tr i -spectrum  will  be  defined  as: 

SxxxK’V'V  =  6  G3(u)l’w2’0j3)  SXxK>  SXX(u)2)  Sxx(w3) 


R  (t  ,t  ,t  )  Expfiw  t  -  i co  t  -  i co  t  )dt  dt  dt 
yxxx'  4’ 5’ 6  14  25  36  456 


f 

* 

1 J 

J 

(59) 


Thus,  if  the  cross-tri-spectrum  can  be  evaluated, an  estimate  for  the 
cubic  frequency  response  function  becomes: 


,  \  Syxxx  ,a)3  ^ 

vwv  =  ro^rrT-CTTTCT 


(60) 


xx'"l'  XX '"2'  XX  3 


This  estimator  is  clearly  analogous  to  that  for  the  quadratic  response, 
Equation  kS.  As  in  the  case  of  the  cross-bi-spectrum,  the  cross-tri¬ 
spectrum  cannot  be  appreciable  unless  the  excitation  spectrum  is 
appreciable  at  all  three  frequencies. 
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By  analogy  with  the  development  of  Equations  50  through  52,  a 
special  fourth  order  correlation  function  may  be  formed  by  replacing 
the  response,  Y(t),  in  Equation  54  with  the  cube  of  excitation: 

X ( t  -  tj)  X ( t  -  t2)  X ( t  -  t3) 


*x3xxx^tl,t2’t3)  =  E 


x3(t)  -  x^TtT 


-  X ( t  -  t  )  x(t-t)  X(t-t  ) 

1  2  3 


-  X(t  -  t2)  X(t  -  tx)  X ( t  -  t3) 

-  X ( t  -  t3)  X(t  -  t?)  X(t  -  t\) 


(61) 


Taking  the  triple  Fourier 

Sx3xxx(w1’W2^3)  " 


transform  and  reducing: 
j||  RX3XXX(tl’t2*t3)  * 

Exp  (-iwjtj  -  i uj2 t o  "  iu3t3)dt1dt2dt3 


=  6  S  (w  ) 
xx  1 


V“2> 


s  (»,) 

xx  3 


(62) 


Accordingly,  an  alternate  identification  form  for  the  cubic  frequency 
response  function  becomes: 


G,  (oj 


■  B,  ,U) 

l’  2’  3 


)  « 


S  (u).  ,u)0  ) 

yxxx  1  2  3 

^X3XXX^1,U)2>w3^ 


(63) 


Thus,  as  in  the  quadratic  case,  the  cubic  frequency  response  function 
may  in  principle  be  estimated  by  the  ratio  of  cross-tr i -spectra . 

It  may  be  noted  that  there  is  a  very  significant  difference  in 

form  between  the  fourth  order  correlation,  R  (t.,tc,t,),  Equation  54, 

’  yxxx  5  6 

and  the  second  and  third  order  correlations,  R  (t)  and  R  (t,  ,  t,), 

*  yx  yxx  5' 

Equations  39  and  45.  This  is  the  presence  of  the  terms  of  the  form: 

Rxx(tj  -  ‘k>  V(t)  X(t  -  •<> 

which  is  expanded  in  Equation  56.  When  the  triple  Fourier  transform  of 


i 
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Equation  56  is  taken  and  a  reduction  by  means  of  the  Parseval  formula, 
Appendix  B  is  carried  out  the  results  are  of  the  form: 


r..  S  (w.)  S  (u.) 
jk  xx  j  xx  l 


MV  +  h 


G3(v,-v,w  )  Sxx(v)dv 


where  F.,  =1  if  u.  =  —to, 
jk  J  k 

=  0  otherwise. 

Thus  the  triple  transform  of  the  fourth  order  correlation  analagous  to 
the  second  and  third  order  correlations,  Equation  53,  contains  delta 
functions  and  the  subtractive  terms  in  the  correlation  function  expressed 
in  Equation  5^  are  for  the  purpose  of  canceling  them  out.  With  respect 
to  the  tr i -frequency  space  in  which  the  cross-tr i -spectrum  is  defined 
the  delta  functions  lie  in  the  planes: 


CO  =  “01 

1  3 

W2  =  "W3 

and  if  they  are  not  removed  somehow  cross-tr i -spectra  1  estimates  of 
the  form  Syxxx(v,-v,w)  will  be  distorted,  perhaps  very  badly. 

If  the  qualitative  results  from  the  simulation  are  applicable 
in  general,  the  surface  defined  by  tri  frequencies  of  the  form  (v,-v,w) 
would  involve  the  portion  of  the  cubic  frequency  response  function 
of  most  practical  importance.  Accordingly  this  portion  of  the 
function  must  be  dealt  wi th, whatever  unpleasant  turns  the  delta  functions 
make  in  the  numerical  analysis. 

The  foregoing  was  about  as  far  as  identification  theory  could 
be  pushed  in  the  present  work.  Several  similar  avenues  of  approach 
(other  forms  of  correlation  functions)  were  attempted,  but  with  no 
better  success,  and  in  general  the  same  form  of  result.  Most  impor¬ 
tantly,  no  relatively  clean  approach  to  the  identification  of  the  linear 
response  function  could  be  developed.  The  conclusion  with  respect  to 
making  further  progress  was  to  postulate  the  following  identification 
strategy: 
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a)  Cross-b i -spectra  1  anlaysis  for  the  quadratic  frequency 
response  function,  Equation  bS  or  52. 

b)  Cross-tr i -spectra  1  analysis  for  the  cubic  frequency 
response  function, especial ly  the  values  of 
G3(v,-v,oj),  Equation  60  or  63. 

c)  Carry  out  a  standard  cross-spectral  analysis  and  make  an 
estimate  of  G^u)  with  a  transposed  form  of  Equation  b] : 

oo 

SyX(w)  "  n  |  G3(“’V’-V)  Sxx(v)dV 

G>)  = - ^ ^ -  (6b) 

xx  (w) 

where  the  estimates  of  G  3(v ,  -v  ,cj)  from  step  b  are  used 
to  "correct"  the  cross  spectrum. 
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CROSS-TR I -SPECTRAL  ESTIMATION 

With  the  tentative  strategy  just  noted  the  only  non-trival  new 
computational  problem  was  the  cross- tr i -spec t rum,  Equation  59.  The 
approach  followed  in  the  development  of  procedures  for  the  cross-bi¬ 
spectrum  in  References  1 k  and  15  followed  the  approach  in  Reference  12 
which  was  in  turn  analagous  to  the  Tukey  autocorrelation  approach  to 
scalar  and  cross  spectral  estimation.  The  Fast  Fourier  Transform 
approach  of  Reference  13  was  by-passed  then  because  it  was  initially 
desired  only  to  probe  specialized  parts  of  the  cross-b i -spectrum ,  and 
because  the  writer  better  understood  the  former  approach.  In  retro¬ 
spect,  the  decision  was  probably  unwise.  In  the  event,  it  was  desir¬ 
able  to  compute  the  entire  cross-bi -spectrum,  an  operation  which  was 
computationally  inefficient  with  the  methods  developed  in  Reference  14. 

However,  some  thought  was  given  to  an  adaptation  of  the  autocorrelation 
approach  to  the  cross-tri -spectrum.  The  attraction  of  this  approach  for  the 

cross-bi-spectrum  had  been  that  a  simple  transformation  of  frequency 
variables  realigned  the  computation  along  and  normal  to  the  lines  of 
symmetry  of  the  cross-bi -spectrum.  The  new  frequency  co-ordinate 
system  was  orthogonal  and  the  result  was  a  third  moment  function  which 
was  easier  to  organize  than  that  indicated  in  Equation  A5.  Accordingly, 
some  effort  was  given  over  to  a  search  for  a  frequency  transformation 
for  tr i -frequency  space  such  that  the  new  frequency  axes  would  lie  in 
the  planes  of  symmetry  of  the  cross-tri -spectrum.  A  number  of  simple 
transformations  were  found,  but  none  were  orthogonal,  and  thus  no 
approach  was  found  to  the  minimization  of  the  computational  problems 
that  Equations  5^  and  59  imply. 

There  was  thus  reason  to  consider  the  methods  of  Reference  27 
more  closely  than  had  been  the  case  previously.  This  reference  was  the 
only  one  known  which  considers  spectra  of  higher  order  than  the  bi¬ 
spectrum,  and  in  which  some  tri-spectral  estimates  were  presented. 

27.  Brillinger,  D.R.  and  Rosenblatt,  M.,  Two  Papers:  "Asympotic  Theory 

of  Estimates  of  kth  Order  Spectra,"  and  "Computation  and  Interpreta¬ 
tion  of  k^  Order  Spectra,"  Proceedings  of  an  Advanced  Seminar  on 
Spectral  Analysis  of  Time  Series,  Edited  by  B.  Harris,  October  1566, 
John  Wiley  6  Sons,  New  York. 
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The  description  in  Reference  28  of  a  more  recent  application  of  the 
methods  to  bi-spectra]  computations  illuminated  some  parts  of  the 
procedure  which  seemed  unclear  in  the  original  reference.  No  explicit 
mention  of  cross-bi  or  cross- t r i -spec t ra  is  made  in  Reference  27. 

However  the  basic  assumptions  involve  a  real,  stationary  vector  process 
(X  (t);  a  =  1,...)  in  which  all  moments  are  assumed  to  exist.  In 

3 

terms  of  the  present  notation  it  is  thus  possible  to  put: 

X, (t)  =  Y(t) 

X:(t)  =  X  ( t ) 

so  that  the  methods  of  Reference  27  appeared  to  apply  to  the  present 
problem,  with  a  minor  redefinition  of  terminology. 

The  basic  idea  of  the  higher  order  spectral  estimates  of  Ref¬ 
erence  27  is  that  a  order  spectral  estimate  is  a  weighted  smoothing 
of  a  kth  order  periodgram  over  a  (k-l)th  order  frequency  space.  In 
terms  of  Reference  27: 

Cross  and  scalar  spectra  are  second  order  (k=2) 

Cross-b i -spectra  are  third  order  (k=3) 

Cross-tri-spectra  are  fourth  order  ( k=^+ ) 

The  adaptation  made  of  the  methods  of  Reference  27  for  present 
purposes  will  be  outlined  in  the  order  in  which  the  computations  would 
be  carried  out.  First  it  is  assumed  that  digital  samples  (N-l)it  seconds 
in  length  are  available  for  the  excitation,  X(t),  and  the  nonlinear 
response,  Y(t).  It  is  common  in  experimental  work  that  the  actual 
level  of  zero  excitation  or  response  is  not  known  to  high  precision. 

Thus  correction  of  time  series  to  zero  sample  mean  (and  removal  of 
evident  trends)  is  a  normal  first  step  in  spectrum  analyses  and  was 
assumed  also  to  be  important  in  higher  order  analyses.  The  response, 
Y(t),  in  all  the  correlations  of  the  last  section  appears  in  the  form 
of  an  implied  correction  to  zero  sample  mean.  Accordingly,  the  first 
step  in  the  process  is  to  remove  trends  and  correct  the  response  to 


28.  Lii,  K.S.,  Rosenblatt,  M.,  and  Van  Atta,  C.,  "Bi-spectral  Measure¬ 
ments  in  Turbulence,"  Journal  of  Fluid  Mechanics,  Vol  .  77, 

Part  I,  1976. 
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zero  sample  mean  so  that  the  resulting  time  series  may  be  defined: 


Y(t)  -  Y(t)  -  Y 1  ( n) ,  n  =  0  ...  N  -  1 

Similarly,  correcting  the  excitation  to  zero  sample  mean  defines  the 
ser i es : 

X(t)-X'(n),  n=0  ...  N-  1 

There  are  two  other  simple  functions  of  excitation  noted  in  the  last 
section,  the  squared  excitation  less  its  expectation  and  the  cubed 
excitation  less  its  expectation.  It  appeared  reasonable  in  these 
cases  to  square  and  cube  X'(n)  for  n  =  0  ...  (N  -  1)  and  to  correct 
the  resulting  series  to  zero  sample  mean  so  that  the  time  series 
corresponding  to  the  squared  and  cubed  functions  of  excitation  will 
be  defined: 

X  ‘ '  ( t )  -  Xr(t)  -  W1 (n),  n  =  0  ...  N  -  1 

X 3  ( t)  -  X-(  t)  -  Z' (n) ,  n  =  0  ...  N  -  1 

It  may  be  noted  that  the  four  operations  just  described  are  exactly 
those  which  would  first  be  carried  out  in  an  implementation  of  an 
autocorrelation  approach  to  cross,  cross-bi,  and  cross- tr i -spectra  1 
ana  lyses . 

One  of  the  subsequent  steps  in  the  procedure  is  to  perform  the 
direct  Fast  Fourier  Transform  on  each  of  the  time  series  defined 
above.  "Tapering"  the  data  to  improve  the  shape  of  the  FFT  spectral 
window  and  to  prevent  "leakage"  from  high  to  low  frequencies  and 
vice  versa  is  fairly  conventional  in  scalar  spectrum  analysis  using 
the  FFT,  is  recommended  in  Reference  27,  and  was  assumed  to  be 
worthwhile  in  the  present  instance.  Thus  the  next  step  in  the  procedure 
was  to  taper  each  of  the  four  time  series  just  described  with  the 
Tukey  10%  cosine  taper  as  described  in  Reference  24.  The  effect,  on 
average,  of  applying  this  taper  function  is  a  "loss"  of  12  1/2%  of 
variance.  To  make  an  approximate  correction  for  this  effect  the 
tapered  time  series  were  finally  multiplied  by  1.069.  At  this  stage 
in  the  procedure  the  four  time  series  may  be  respresented  as  follows: 
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V  (t)  -  Y  ( t )  *  Y(n) ,  n  =  0,  1  ...  N  -  1 

X( t)  —  X ( n) ,  n=0,  1  ...  N-  1 

X  ( t )  -  x- (t)  -  W ( n) ,  n  =  0 ,  1  ...  N  -  1 


X  '  (  t)  -  X;  ( t)  -  Z(n) ,  n  =  0,  1  ...  N  -  1 

where  each  tine  series  is  corrected  to  zero  sampie  mean,  tapered,  and 
compensated  tor  loss  of  variance. 

The  next  step  in  the  procedure  is  to  perfo-n  the  direct  Fast 
Fourier  Transform  on  each  series.  The  FFT  algorithm  used  in  the 
present  instance  performs  the  following  operations: 

,  n:1 

Y(k)  =  —  Y(n)  Exp  (-i2'nk/N) 

N  n=0 


(65) 


,  N;' 

X  ( k)  =  -  X ( n)  Exp  ( - i 2-nk/N) 

N  n=0 


(66) 


1  N-1 

W(k)  =  tt  ^  W(n)  Exp  (-i2nnk/N) 
n=0 


(67) 


i  n:1 

Z(k)  =  TT  '  Z(n)  Exp  (-i2nnk/N) 
N  n=0 


(67) 


where  the  transforms  are  defined  for: 


,N 


-(j  -  D  <  k  <  N/2 


circular  frequency  is  related  to  the  index  k  by: 

NAt 


and  the  transform  for  negative  k  is  the  complex  conjugate  of  that  for 
pos i t i ve  va 1 ues : 

Y(-k)  =  Y* ( k ) 

X(-k)  =  X*(k) 

W(-k)  =  WMk) 

Z(-k)  =  Z*(k) 
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Apart  from  the  sequence  of  the  compensat'cn  for  variance  loss,  the  oper¬ 
ation  thus  far  is  exac  'y  that  which  would  be  carried  out  in  deriving 
scalar  spectra  of  X(t),  Y(t),  a' (t)  and  X'(t)  by  the  FFT  method. 

It  is  necessary  next  to  define  the  periodograms  required  by 
the  methods  of  Reference  27  in  present  notation.  Starting  with  the 
second  order  case,  an  obvious  estimate  of  the  cross  correlation  function, 
Equation  39.  for  an  N  point  sample  may  be  written: 

R  (P)  =  i  7  Y(n)  X (n  +  p)  (69) 

yx  N  „-n 

n  L/  's 

where  the  domain  of  the  summation  is: 

0  <  n  <  N  -  1 
0<n-p<N  -  1 

The  second  order  periodogram  is  the  discrete  Fourier  transform  of 
Equation  69: 

N-1 

P  (j)  =  At  l  R  (p)  Exp(-i2irpj/N) 

Y  p—  N+1  YX 

N+1  (.  ) 

=  NAt  V  ^  —  Y(n)  Exp(+i 2nnj/N)  • 

n=02  p=-N+1 

] 

X(n+p)  Exp{*  i  2nj  (n+p)/N)  (70) 

=  NAt  Y(-j)  X(j) 

The  result,  Equation  70,  is,  apart  from  notation,  the  "raw"  cross 
spectrum  of  Reference  2b. 

The  third  order  periodogram  may  be  approached  similarly  by 
defining  an  estimate  of  the  third  order  correlation  function  of 
Equation  ^5  as: 

R  (P.q)  *  tr  [  Y(n)  X(n-p)  X(n-q)  (71) 

"  n=03 
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where  the  domain  D_ i s  defined: 

f  1 

0  <  n  <  N  -  1 

0<n-p<N-1 

0<n-q<N-1 

( 

The  corresponding  third  order  periodogram  is  the  double  Fourier  transform 
of  Equat  ion  71  •' 

N-1  N-1 

P  (J.k)  =  (At)2  I  I  R  (p,q)  Exp  (  -  i  2rpj  /N-  i  2’iqk/N) 

YXX  p=-N+1  q=-N+1  yXX 

=  (NAt)2  Y(j  +  k)  X(-j)  X(-k)  (72) 

where  the  derivation  is  similar  to  that  of  the  second  order  case.  The 
third  order  periodogram  corresponding  to  the  squared  excitation  correla¬ 
tion,  Equation  50,  may  be  written  by  analogy: 

Pwxx(j’k)  =  (N't)2  W(j  +  k)  X(-j)  X(-k)  (73) 

Next,  the  fourth  order  periodogram  is  approached  by  defining  a 
fourth  order  correlation  estimate  analogous  to  the  first  term  of 
Equation  5^: 
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The  corresponding  fourth  order  periodogram  is  the  triple  Fourier  transform 
of  Equation  lk: 


P 

yxxx 


( j  ,  k,x)  =  (At)  ' 


N+1 

p=-N+ 1 


N+1 

q=-N+) 


N+1 

r=-N+ 1 


R 


yxxx 


(p,q, 


r)  • 


Exp  Q  i  2*pj  /N-  i  2-qk/N-  i  2”  r  .  /fj] 

•  =  (NAt)3  Y(j  +  k  +  i)  X(-j)  X(-k)  X(-0  (75) 

and  the  periodogram  corresponding  to  the  first  term  of  the  cubed 
excitation  correlation.  Equation  61,  may  be  written: 

» 

Pzxxx(j,k,Z)  =  (NAt)3  Z(j  +  k  +  0  X(-j)  X(-k)  X(-«)  (76) 

It  may  be  noted  that  the  frequency  indices  in  the  final  expressions 
for  the  per iodog rams  (Equations  70,  72,  73,  75,  76)  sum  to  zero.  That 
is,  in  Equation  75  for  example, 

(j  +  k  +  l)  +  (-j)  +  (-k)  +  {-l)  =  0 

This  is  a  property  required  by  the  definition  of  the  kt*1  order  periodogram 
of  Reference  27. 

It  was  shown  in  Reference  27,  for  the  strictly  stationary  processes 
I  presently  assumed,  that  the  expected  value  (or  temporal  mcon)  of  the  k*^ 

order  periodogram  equals  the  k*^  order  spectrum,  but  with  an  important 
provision.  The  provision  has  to  do  with  the  presenceof  delta  functions 
discussed  in  the  last  section.  Essentially,  as  sample  length  approaches 

•  infinity,  the  periodogram  approaches  the  spectrum  except  in  those 
particular  regions  of  frequency  space  where  delta  functions  are  expected. 
However  it  was  suggested  that  reasonable  estimates  may  be  made  by  a 
multi-dimensional  frequency  smoothing  operation  on  the  periodograms  if 

>  the  problem  regions  are  avoided. 


I 
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As  in  lower  order  spectrum  analyses  from  finite  samples  it  is 
possible  only  to  estimate  higher  order  spectral  averages  of  the  form: 


S(u)j,u)0,.#)  —  •  .  W(aJj  “ 


)  ,n:...)d^:di:2. 


where  the  weighting  function  W(a  peaks  at  frequency  vector 

(0,0,...),  falls  off  to  zero  in  every  direction,  and  has  the  property: 


1  =jj..j  W(a1,o0...)  da i da„ . . . 

No  specific  guidance  is  given  in  Reference  27  as  to  what  the  discrete 
analog  of  the  weighting  function  ought  to  be  for  the  higher  order  case. 
By  analogy  with  conventional  frequency  weighting  (smoothing)  techniques 
for  scalar  and  cross  spectra,  Reference  24,  a  reasonable  solution  was 
anticipated  if  W(a,,a„...)  was  zero  everywhere  except  within  some 
reg ion : 


-e  <  a1  <  e 
-c  <  a „  <  e 


and  a  constant  otherwise  (essentially  a  multiple  dimension  block 
average).  This  was  found  to  be  the  case  in  the  analysis  of  Ref¬ 
erence  28,  and  thus  encouraged,  the  block  average  technique  was 
adopted  for  the  present. 

Wi th  this  decision,  the  estimating  forms  for  the  various  cross 
spectra  of  present  interest  could  be  translated  from  the  basic 
development  of  Reference  27.  The  cross-spectral  estimate  correspond¬ 
ing  to  Equation  38  becomes: 


kl  A  ..  j+m  _  _ 

s  (j)  =  l  Y(-h)  X(h)  4(h) 

h-j  “m 

where,  as  before,  circular  frequency  is  defined  as: 

-  2jLl 
"  NAt 
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"m"  is  the  block  size  parameter  corresponding  to  e,  "H"  equals  the  number 
of  terms  summed,  and 

e  j  ( h )  =  0  ifh  =  0 
=  1  otherwise 

The  switch  function  $  (h)  suppresses  a  delta  function  at  zero  frequency. 
Otherwise,  the  estimate,  Equation  77,  is  exactly  the  conventional 
FFT  frequency  smoothed  cross  spectral  estimate. 

Cross-b i -spectra  1  estimates  corresponding  co  Equations  kb  and  k8 

become: 


> 


9 


» 


I 


» 


„  /...  »2  J+m  k+m  _ 

s  (j,k)  =  -  I  I  Y(f  +  h)  X(-h)  X(-f)  *  (h,f)  (78) 

yXX  H  f-j-m  h=k-m  2 

-  j4™  l<+m  _  _ 

Sx2xx(j,k)  =  ~H  '  I  l  W(f  +  h)  X(-h)  X(-f)  4*2  (h,f)  (79) 

A  f -j -m  h=k-m 

The  switch  function  2 ( h , f )  is  defined: 


<J>2  (h,  f )  =0  if:  h  =  0 
or:  f  =  0 
or:  h  +  f  =  0 
<f>2(h,f)  =  1  otherwise 

In  this  case  it  may  be  noted  that  the  development  of  Reference  27  predicts 
delta  functions  along  the  line  j  =  -k  in  the  bi-frequency  plane,  and 
the  switch  function  suppresses  the  periodogram  estimates  along  this  line. 
This  is  in  contrast  to  the  autocorrelation  based  methods  of  Reference  )k 
and  15  with  which  no  such  problem  appeared  to  exist. 


Finally,  cross-tr i -spectra  1  estimates  corresponding  to  Equations  59 
and  62  become: 


„  u  j+m  k+m  £+m  _ 

Jj.M)  =  ~ J  ~  l  l  I  V  (e+f+h)  X(-e)X(-f)X(-h)<J>(e,f  ,h) 

yxx  eaj-m  f=k-m  h=t-m 

(80) 
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s  , 

X  XXX 


(j ,k, i) 


(NAt)~ 

H 


j+m  k+tn  i+m  _ 

l  T  Z(e+f+h)X(-e)X(-f)X(-h) 

e=j-m  f=k-m  h=« -m 


The  switch  function  ;  (e,f,h)  is  defined: 


^ 3 (e, f , h) 


V-,  (e,f  ,h) 


=  0  if: 
or : 
or : 
or : 
or : 
or : 
or: 


e  =  0 
f  =  0 
h  =  0 
e+f+h  *  0 
e+f  =  0 
e+h  =  0 
=  0 


f+h 

1  otherwise 


As  would  be  anticipated  from  the  discussion  of  the  last  section, 
the  switch  function  suppresses  the  delta  functions  in  the  periodogram 
which  lie  exactly  in  the  planes  of  the  cross-tri-spectrum  which  may  be 
of  considerable  practical  importance.  As  previously  noted,  i t  was  hoped 
that  averaging  over  (2m)  adjacent  periodogram  estimates  would  provide 
reasonable  results  in  these  cases. 

As  in  conventional  spectrum  analyses  it  is  expected  that  there 
is  little  point  in  producing  estimates  of  any  of  the  cross  spectra  for 
frequencies  spaced  at  intervals  less  than  Z^(2m  +  1)/N£t.  The  analysis 
parameter  "m"  controls  the  extent  of  frequency  smoothing  and  must  be 
chosen  in  accordance  with  the  problem  at  hand  and  the  frequency  resolu¬ 
tion  of  the  basic  Fast  Fourier  transforms.  So  long  as  the  minimum 
frequency  increments  correspond  to  (2m  +  l)  there  appears  no  disadvant¬ 
age  in  basic  programming  suitable  for  probing  the  cross-bi  and  cross-tri 
spectra  as  opposed  to  a  complete  computation  for  all  combinations  of 
frequency. 

The  primary  reason  for  the  frequency  smoothing  operations  is 
to  reduce  the  variance  of  the  estimates.  There  appears  no  reason  in 
References  27  and  28  why  simple  ensemble  averaging  of  estimates  should 
not  also  be  carried  out  to  the  extent  that  multiple  independent  time 
domain  samples  are  available. 

8A 


; ,  ( e ,  f ,  h ) 
(81) 
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CROSS-B! -SPECTRAL  IDENTIFICATION  OF  THE  QUADRATIC 
FREQUENCY  RESPONSE  FUNCTION 

The  last  objective  of  the  present  work  was  to  attempt  the  appli¬ 
cation  of  the  approaches  described  in  the  last  two  sections  to  the 
simulated  random  data.  As  noted  previously,  the  simulated  response  data 
had  been  produced  in  such  a  way  that  response  could  easily  be  computed 
for  various  excitation  levels.  Accordingly,  the  first  decision  to  be 
made  was  what  excitation  level  to  employ.  It  appeared  from  the  analysis 
of  scalar  spectra  that  the  response  of  the  simulated  system  to  excitation 
twice  the  nominal  was  in  line  with  nonlinear  response  which  might  be 
observed  in  towing  tank  experiments,  and  this  level  (o  =  0.5)  was 
selected.  (The  corresponding  scalar  spectra  are  shown  in  Figure  lA.) 

For  convenience,  the  initial  basic  steps  noted  in  the  previous 
section  were  carried  out  on  each  of  the  eleven  samples  of  excitation 
and  response  which  were  available,  and  the  results  were  stored  for  sub¬ 
sequent  use.  In  particular,  the  processing  steps  carried  out  were: 

1.  Compose  the  samples  to  be  analyzed  by  multiplying  the 
nominal  excitation  series  and  response  component  Y1 (n) 
by  two,  Y2(n)  by  A,  and  Y3(n)  by  8.  Summing  the 
response  components  in  accordance  with  Equation  lA, 
there  resulted  11  sets  of  time  history  of  response, 

Y(n)  and  excitation  X(n). 

2.  Correct  both  series  to  zero  sample  mean,  and  derive 
W'(n)  and  Z * ( n )  as  previously  described. 

3.  Taper  each  series  and  correct  for  variance  loss. 

A.  Perform  the  direct  Fast  Fourier  Transform  on  each 

series,  so  that  X(k),  Y(k),  W(k)  and  Z(k)  were  avail¬ 
able  as  in  Equations  65  through  67. 

The  first  step  of  the  analysis  strategy  previously  noted  is  to 
perform  cross-b i -spectra  1  analvses  to  define  quadratic  response.  Since 
this  type  of  analysis  has  been  done  previously,  there  was  no  particular 
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question  of  feasibility  involved.  However,  in  the  present  instance  it 
was  though  worthwhile  to  perform  the  cross-bi-spectral  analysis  with  the 
FFT  based  method  to  gain  some  experience  and  to  see  if  the  switch  func¬ 
tion  shown  after  Equation  79  would  adequately  cope  with  the  delta  functions 
expected  along  the  line  y_,  =  -*j  in  the  bi-frequency  plane.  To  this 
end  Equations  78  and  79  were  programmed  together  in  order  to  save  some 
complex  multiplications,  and  applied  to  the  11  sets  of  transform  data 
previously  described.  The  block  averaging  parameter  "m"  was  set  at  b; 
that  is,  81  adjacent  values  of  the  third  order  periodogram  were  averaged 
to  form  the  cross-bi-spectral  estimate  at  frequencies  =  j£w,  u)„  =  k£u 
where 

Aoj  =  =  0-^22  rad/sec 

NAt 


The  result  was  a  set  of  cross-b i -spect ra 1  estimates  for  each  of  the  11 
samples  of  data,  and  the  final  estimates  were  made  by  ensemble  smoothing; 
that  is,  by  averaging  the  11  sets  of  cross-bi-spectral  estimates. 


Qualitatively,  the  resulting  cross-bi-spectral  estimates  appeared 
reasonable,  and  it  appeared  that  the  most  sensitive  way  to  check  the 
result  would  be  to  form  an  estimate  of  G  (wj.oo,,)  and  compare  the  results 
with  the  known  theoretical  quadratic  frequency  response  function.  To 
this  end  the  quadratic  frequency  response  function  was  estimated  by: 


_  S  (u),  ,  U) ^  ) 

„  /  ^  yxx  l’ 

-  t - 


S  •>  (u.  ,a,) 

X“XX  1  2 


where  the  cross-b i -spectra  were  the  final  ensemble  smoothed  values. 

Previous  experience  with  this  type  of  analysis  has  suggested  that  such 
estimates  cannot  be  good  outside  the  range  of  bi-frequency  where  the 
product  S  (w  )  S  (w.)  is  appreciable.  Inspection  of  S  -  (w.,w,) 

XX  1  XX  X  XX  l  * 

suggested  that  only  in  the  range: 

5  <  |cjj  |  <  10 

5  <  M  <  10 

would  this  be  reasonably  satisfied,  and  accordingly  attention  was  paid 
only  to  this  range  of  bi-frequency.  For  comparison,  theoretical  values 
of  G2(u)j,u)2)  were  computed  by  means  of  Equations  11  and  13  for  the  discrete 
bi-frequencies  of  the  estimates  from  the  simulated  samples. 
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The  resulting  comparisons  are  shown  in  Figures  18  through  20. 

Figure  18  involves  the  modulus  of  Grfu^.w.,),  Figure  19  the  real  part 
and  Figure  20  the  imaginary  part.  These  figures  are  essentially 
plotted  tables.  In  the  region  of  the  bi-frequency  plane  where  both 
frequencies  are  positive  and  the  frequencies  are  between  5  and  10  rad/sec 
the  procedure  just  described  results  in  a  matrix  of  121  estimates. 

Because  of  the  symmetry  of  the  function  it  suffices  to  show  only  the 
65  estimates  where  Oj  >  u)2 ,  and  these  are  plotted  at  their  location 
in  the  plane  in  vertical  lettering.  Just  below  each  of  the  estimates 
the  corresponding  theoretical  value  is  shown  in  slanted  lettering. 

Similarly  65  unique  estimates  of  G,(uJ1,aj?)  result  for  the  case  that 
oi  and  ui2  differ  in  sign  and  these  are  shown  in  the  lower  part  of  the 
figures  along  with  the  corresponding  theoretical  values. 

In  general,  the  correspondence  between  theoretical  and  estimated 
G-,(wj,i»>2)  shown  in  Figures  18  through  20  seems  quite  good,  in  fact, 
better  than  might  have  been  expected  on  the  basis  of  past  experience 
with  samples  of  about  the  same  size  and  the  cross  correlation  approach. 
References  14  and  15-  Signs  and  magnitudes  of  real  and  imaginary  parts 
are  in  quite  reasonable  agreement  except  for  the  values  shown  at 
w1,oi2  =  (9.7,  -9.7).  In  this  instance  Sx2xx(wi>w2^  ‘s  verV  much  smaller 
than  adjacent  estimates,  and  it  was  suspected  that  the  bad  result  was 
largely  a  matter  of  division  by  noise.  It  is  along  the  line  =  -w 
in  the  lower  part  of  the  figures  where  wild  disagreement  due  to  the 
presence  of  delta  functions  would  be  expected.  With  the  exception  just 
noted,  percentage  differences  between  estimate  and  theory  along  this 
line  do  not  seem  greatly  different  from  those  off  the  line.  The  net  result 
of  the  cross-bi-spectral  exercise  was  to  suggest  that  cross-tri-spectral 
estimation  should  next  be  attempted  since  the  switch  function  in 
Equations  78  and  79  appeared  to  be  coping  with  the  anticipated  problem 
area. 
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Estimates:  Vertical  Lettering 
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TRIAL  CROSS-TRI -SPECTRAL  IDENTIFICATIONS 


Since  the  main  objective  was  to  identify  the  cubic  frequency  response 

function  from  simulated  data,  end  the  theoretical  response  function  was 

all  that  was  known,  the  trial  analyses  followed  the  general  approach  of 

the  last  section.  That  is;  the  estimates  defined  by  Equations  80  and  8l 

were  programmed  together  so  as  to  allow  computation  of  the  cross-tri- 

spectra!  estimates  S  (w,  )  and  S  ?  ,w,)  ,  and  then  the 

yxxx  1-3  xxxx  *  - 

cubic  frequency  response  function  was  estimated  by: 


r  i  \  _  yxxx  - 


S  3  (a.1,  ) 

x  xxx v  1  2  3 


in  accordance  with  Equation  63.  Trial  computations  of  theoretical  values 
of  G3 ,w0 ,w3)  suggested  that  the  frequency  resolution  of  the  analysis 
should  not  be  coarser  than  that  of  the  quadratic  analysis  of  the  last 
section  and  the  same  block  averaging  parameters  were  picked.  The  block 
averaging  parameter  "m"  was  thus  set  at  4,  and  accordingly  729  adjacent 
estimates  of  the  fourth  order  per iodog rams  were  to  be  averaged  in  order 
to  form  the  cross- tr i -spect ra 1  estimates  at  tr i -f requenc i es  which  would 
be  multiples  of  Aw  =  0.442.  The  same  basic  transform  data  (X(k),  Z(k) 
and  V(k))  produced  for  the  cross-b i -spectra  1  analysis  was  required  for 
the  cross-tri-spectra! .  Thus  the  starting  point  for  the  present  analysis 
was  the  eleven  sets  of  transforms  previously  used  which  correspond  to 
the  simulations  for  a =  0.5  (spectra  in  Figure  14).  As  in  the  cross-bi- 
spectral  case, ensembl e  smoothing  was  performed  over  the  cross-tri-spectra  1 
estimates  formed  from  each  of  the  11  simulated  data  records. 


It  is  clear  from  the  basic  identification  theory  that  cross-tri 
spectral  estimates  cannot  be  appreciable  outside  the  range  of  tri¬ 
frequency  where  the  product: 


S 

xx 


(w  )  S 

1  xx 


)  S  UJ 

XX  - 


is  appreciable.  Estimates  of  G  (uu^  ,w3 )  outside  this  range  are  apt  to 
be  seriously  corrupted  by  noise,  Some  preliminary  computations  of 


> 


v— yv«f* 
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S  3  (u)  ,jj  ,d)  )  and  Inspection  of  the  excitation  spectrum  of  Figure 

X  XXX  1  Z  3 

suggested  that  only  in  the  range 

5  <  l&jjl  <8  radians/second 

5  <  |a.2|  <  8 

5  <  |u3|  <  8 

would  there  be  hope  of  reasonable  results. 

Since  an  analysis  over  even  a  restricted  region  of  tri-frequency 
space  would  require  appreciable  computation  time,  it  was  determined  to 
start  slowly  and  first  examine  some  simple  cases.  The  simplest  case  is 
just  to  estimate  (w,w,w) ,  the  third  harmonic  response  function.  In  this 
case  delta  function  problems  are  not  of  concern.  Accordingly,  estimates 
of  cross-tr i-spectra  along  the  line  *  u>2  *  *  w)  were  formed  from 

each  of  the  eleven  samples,  averaged,  and  from  these  results  estimates 
of  the  cubic  response  function  G3 (w.to.w)  were  made.  The  results  are 
shown  as  circled  points  in  Figure  21.  In  the  figure  the  real  and 
imaginary  parts  as  well  as  the  modulus  of  G3(w,w,w)  are  plotted  for  the 
restricted  range  of  frequency  just  noted.  The  corresponding  theoretical 
values  of  G3(u),w,u))  are  shown  as  dashed  lines.  It  will  be  noted  that 
arrows  have  been  drawn  to  indicate  the  direction  of  the  estimates  made 
for  a  frequency  just  below  5  rps  (the  positions  of  these  points  are  off 
the  plotti ng  field) . 

It  may  be  noted  from  Figure  21  that  the  results  of  this  first 
cross-tri-spectral  identification  (circled  points)  were  at  least  recog- 
niszble  as  an  estimate  of  the  theoretical  values  of  the  function. 

Clearly,  significant  deviations  were  present.  Those  at  the  ends  of  the 
freauency  range  are  typical  of  what  can  happen  when  the  excitation  cross- 
tri-spectrum  gets  small.  However  the  deviations  in  mid  range  gave  rise 
to  the  thought  that  perhaps  the  sample  size  (11  records)  was  too  small. 
There  is  no  theory  of  sampling  variability  available  for  cross-tri-spectral 
estimation  of  cubic  response  functions.  (There  is  none  for  cross-bi- 
spectral  identification  either.) 


92 


u,  Radlans/Sacond 


FIGURE  21  COMPARISON  OF  ESTIMATED  AND 
THEORETICAL  VALUES  OF  THE 
CUBIC  FREQUENCY  RESPONSE  FUNCTION 
ALONG  THE  LINE  w,  ■  ti.  ■  ui,  *  co. 


TP-2275 


An  obvious  next  step  was  to  generate  significantly  more  sample  and 
repeat  the  identification.  It  was  relatively  convenient  to  modify  the 
programs  described  in  conjunction  with  the  basic  simulation  so  that  an 
additional  20  statistically  independent  time  domain  sample  records  of 
excitation  could  be  generated.  Computation  of  the  time  domain  responses 
for  each  of  these  was  carried  out  as  before,  with  the  result  that  the 
simulated  data  set  was  increased  from  11  to  31  records — roughtly  tripling 
the  available  sample.  The  cross-tr i-spectral  analysis  was  repeated  with 
ensemble  smoothing  over  the  31  records,  and  the  results  are  shown  in 
Figure  21  as  triangles. 

The  results  in  Figure  21  for  31  records  are  clearly  an  improvement 
over  the  first  analysis.  Agreement  between  estimate  and  theory  is 
reasonably  good  at  the  high  end  of  the  frequency  range.  At  the  low  end 
of  the  range  the  deviations  between  estimate  and  theory  were  reduced, 
but  not  by  a  great  deal.  The  results  of  varying  sample  size  suggest 
that  the  point  of  diminishing  returns  may  be  somewhere  between  11  and 
31  records  so  that  the  idea  of  generating  even  more  sample  was  discarded. 

It  has  been  noted  that  the  third  harmonic  response  of  the  present  system 
is  quite  small.  In  this  light  the  results  in  Figure  21  might  be  considered 
not  too  bad  for  a  first  attempt. 

The  next  set  of  trial  identifications  were  for  the  cubic  response 
function  along  the  line: 

“l  *  W2  *  “ 

“3  *  "u 

This,  according  to  the  speculations  previously  made  may  generally  be  the 
most  practically  important  part  of  the  cubic  frequency  response  function. 
G3(u,b>,-u)  expresses  cubic  nonlinear  response  at  excitation  frequency, 
and  is  part  of  the  special  case,  G3(w,v,-v),  which  must  be  identified 
if  the  linear  response  function  is  to  be  separated.  Cross-tri-spectral 
identification  was  carried  out  as  before  for  this  case  with  the  results 
shown  in  Figure  22.  This  figure  indicates  the  real  and  imaginary  parts 
as  well  as  the  modulus  of  estimated  and  theoretical  values  of  G^ (u>,w,-ui) . 
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Ensemble  smoothed  estimates  for  11  and  31  records  are  shown  by  circular 
and  triangular  symbols.  As  in  the  previous  trial  the  differences 
between  estimates  for  11  and  31  samples  imply  that  sample  size  is 
acceptable. 

Comparing  the  estimates  and  the  theory  in  Figure  22,  it  is 
clear  that  something  is  wrong.  The  magnitude  of  the  estimates  are  in 
the  ballpark  but  the  signs  of  the  real  and  imaginary  parts  are  consis¬ 
tently  opposite  what  they  should  be. 

An  obvious  diagnostic  step  was  to  see  if  the  presence  of  a  quad¬ 
ratic  component  (Y2(n),  Equation  14)  was  upsetting  things.  To  this  end 
new  time  domain  samples  with  quadratic  component  suppressed  were  made 
up  for  a  few  records  and  the  identification  repeated.  Omitting  the 
quadratic  component  had  little  effect  on  the  estimates  from  the  cross- 
tri-spectral  identification.  This  is  of  course  what  is  expected  from 
the  theory. 

Continuing  the  diagnosis,  new  time  domain  samples  which  contained 
only  the  cubic  component  (Y(3)  of  Equation  14)  were  constructed  and  the 
cross-tr i-spectral  identification  repeated.  The  results  are  shown  by 
rectangular  symbols  in  Figure  22.  (Ensemble  smoothing  over  all  31 
time  domain  samples  was  involved.)  The  estimates  of  G3(w,w,-u>)  for 
the  case  where  only  the  cubic  component  is  involved  in  the  data  are  quite 
reasonable  relative  to  the  theoretical  values. 

It  is  obvious  from  this  evidence  that  there  is  a  problem  with 
the  delta  functions,  since  they  will  be  involved  in  this  case.  The 
theory  of  Reference  27  says  essentially  that  reasonable  results  should 
be  obtained  for  averages  near,  "but  not  too  near"  the  delta  function  tri¬ 
frequencies.  It  was  thought  possible  that  the  influence  of  the  deltas 
was  spread  out  over  a  few  more  periodogram  estimates  than  those  which 
are  exactly  on  the  critical  tr 1-frequencies.  Accordingly,  the  switch 
functions,  Equations  80  and  8l ,  were  modified  so  that  periodogram 
estimates  close  to  but  not  at  the  critical  points  could  also  be  neglected. 
The  modified  program  was  exercised  so  as  to  do  the  identification  with 
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neglect  of  successively  wider  bands  of  periodogram  estimates.  No 
essential  change  in  the  character  of  the  estimates  was  noted  up  to  a 
neglect  band  which  amounted  to  almost  half  the  delta  frequency  of  the 
analysis.  At  this  point  the  scatter  of  estimates  increased  radically 
since  far  fewer  periodogram  estimates  were  being  averaged.  The  con¬ 
clusion  is  that  the  influence  of  the  delta  functions  "spreads  out" 
considerably  from  the  exact  critical  tri-frequencies. 

Clearly,  the  presence  of  the  linear  component  in  the  data  strongly 
influences  the  importance  of  the  delta's  (see  discussion  after  Equation  63). 
It  appeared  that  the  cross-tri -spectral  identification  method  developed 
would  only  yield  reasonable  results  for  G3 (w.w.-w)  when  the  linear  com¬ 
ponent  was  relatively  small;  that  is,  for  excitation  levels  in  which  the 
major  part  of  the  response  was  cubic.  This  is  not  likely  to  happen 
in  practical  ship  seakeeping  problems  so  that  the  method  as  developed 
appears  unsatisfactory  for  G3(a),u),-&>) .  Since  G3(w,w,-w)  is  probably 
the  most  important  part  of  G3(w,v,-v)  the  method  is  almost  certainly 
an  unsatisfactory  approach  to  the  latter  case.  There  seemed  no  point 
in  proceeding  further  until  a  method  which  was  satisfactory  for 
Gg (ui,(ij,-oj)  was  in  hand,  since  indentif ication  of  the  linear  frequency 
response  function  depends  upon  success  with  the  identification  of  cubic 
response  function  of  this  type. 
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CONCLUDING  REMARKS 

The  objectives  of  the  present  work  were  to  explore  the  applicability 
of  the  third  degree  functional  polynomial  model  to  nonlinear  seakeeping 
problems,  and  to  attempt  the  development  of  analysis  approaches  by  which 
third  degree  nonlinearities  in  the  responses  of  ships  in  waves  might  be 
interpreted  and  characterized. 

It  appears  that  the  functional  polynomial  of  third  degree  (which 
contains  linear,  quadratic  and  cubic  terms)  may  be  capable  of  representing 
a  wide  variety  of  the  relatively  weak  nonlinearities  which  can  be  antici¬ 
pated  in  the  response  of  ships  to  waves.  Conceptually,  the  ship  dynamics 
would  be  represented  by  a  series  of  frequency  response  functions--! inear, 
quadratic  and  cubic.  These  functions  can  be  interpreted  in  terms  of  the 
responses  to  the  superposition  of  one,  two  and  three  periodic  excitations. 
Thus  in  principle  the  functions  may  be  identified  by  deterministic 
experiments.  It  appears  that  extension  of  hydrodynamic  theory  to 
include  cubic  non  1  inear i ties  may  be  accomplished  at  the  expense  of 
systematic  development  to  "third  order"  just  as  it  was  found  possible 
to  include  quadratic  nonlinearities  with  a  development  to  "second  order". 

Within  limits,  the  functional  polynomial  model  is  compatible 
with  the  random  excitation  case.  Past  work  has  shown  the  qualitative 
and  quantitative  influence  of  quadratic  nonlinearities  upon  the  response 
to  random  excitation.  A  significant  part  of  the  present  work  was  to 
explore  the  influence  of  cubic  nonlinearities.  This  was  done  by  means  of 
a  simulation  of  a  relatively  arbitrary  nonlinear  system  so  that  the 
indications  are  hardly  general.  However  in  many  respects  the  response 
of  the  simulated  system  reflected  some  types  of  nonlinear  behavior  which 
has  been  seen  in  experiment  so  that  some  reasonable  conjectures  can  be 
made.  The  first  is  that  the  most  significant  aspects  of  cubic  nonlinear 
time  domain  response  to  random  waves  may  be  associated  with  wave  groups, 
just  as  has  been  noted  in  the  quadratic  case.  The  second  is  that  the 
most  important  influence  of  cubic  nonlinearities  on  the  scalar  spectrum 
of  nonlinear  response  appears  to  be  dependent  upon  a  special  portion  of 
the  cubic  frequency  response  function.  This  special  portion  of  the 
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function  involves  third  degree  nonlinear  interactions  between  two 
frequency  components  as  well  as  a  self  interaction  involving  a  single 
frequency. 

It  was  important  in  the  development  of  credibility  of  both  the 
linear  and  the  quadratic  functional  models  as  applied  to  ship  response 
that  means  be  developed  of  identifying  linear  and  quadratic  frequency 
response  functions  from  samples  of  random  excitation  and  response. 
Accordingly,  a  considerable  part  of  the  present  work  was  given  over  to 
the  development  of  theoretical  identification  methods  for  the  cubic 
response  function.  The  theoretical  developments  in  this  direction 
included  the  definition  of  an  entity  named  the  cross-tri-spectrum, 
which  is  simply  related  to  the  cubic  frequency  response  function,  and 
thus  in  principal  to  an  identification  method. 

The  theoretical  development  disclosed  a  complication  with  cubic 
systems  not  present  in  the  linear  plus  quadratic  model.  This  is  that 
the  linear  and  cubic  responses  are  mixed  together  in  such  a  way  that 
the  results  of  a  conventional  cross  spectrum  identification  cannot  be 
simply  related  to  the  linear  response  function.  To  the  extent  that 
it  was  possible  to  progress  the  theory,  it  appears  that  in  analyzing 
linear  plus  cubic  systems  a  portion  of  the  cubic  response  function 
must  be  identified  first.  The  theory  also  disclosed  that  the  cross- 
tri-spectrum  is  not  a  simple  extension  of  cross  and  cross-b?-spectral 
theory.  The  difference  is  related  to  the  problem  with  the  cross 
spectrum.  Essentially,  unless  special  provision  are  made,  a  cross-tri¬ 
spectrum  will  contain  delta  functions  which  are  strongly  influenced 
by  the  linear  components  of  response.  Unfortunately,  the  location 
of  these  delta  functions  corresponds  exactly  to  that  of  the  portion 
of  the  cubic  frequency  response  function  thought  to  be  of  most  practical 
importance. 

The  last  part  of  the  present  work  was  to  attempt  the  cross-tri- 
spectral  identification  of  the  cubic  response  functions  from  simulated 
random  excitation  and  response  data.  There  were  two  avenues  of  approach 
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available  for  cross-tr i-spectral  estimation,  a  correlation  function 
approach  thought  to  be  extremely  expensive,  and  a  Fast  Fourier  Transform 
approach.  The  FFT  approach  was  opted  for  and  both  cross-bi  and  cross-tr 
spectral  estimating  methods  were  developed.  Identifications  of  the 
quadratic  frequency  response  function  were  reasonably  successful  with 
this  approach.  However  when  trial  identifications  of  cubic  frequency 
response  functions  were  carried  out  troubles  surfaced  immediately. 
Essentially,  the  methods  developed  appear  workable  for  those  portions 
of  the  cross-tri-spectrum  where  no  delta  functions  are  expected,  and 
for  the  portions  which  have  delta  functions  when  no  linear  component 
of  response  is  present.  Thus  the  FFT  methods  developed  were  found  un¬ 
satisfactory  for  the  portion  of  the  cubic  frequency  response  function 
of  most  practical  importance,  so  that  there  appeared  little  point  in 
further  pursuit  of  the  FFT  based  cross-tr i-spectra I  analysis  method. 
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,  RECOMMENDATIONS 

Clearly,  the  FFT  cross-tri -spectral  estimating  approach  attempted 
in  the  present  work  appears  a  failure  in  the  case  of  most  practical 
I  importance.  This  does  not  necessarily  mean  that  the  basic  identification 

theory  is  incorrect  or  unusable,  because  the  FFT  approach  was  selected 
largely  upon  economic  grounds.  If  the  third  degree  functional  model  is 
to  achieve  credibility  some  means  of  identification  from  random  wave 
p  data  will  be  required.  The  obvious  next  step  in  the  development  of 

cross-tri-spectral  analysis  is  to  go  back  to  the  theoretical  correlation 
function  approach  and  to  develop  numerical  estimation  methods  on  that 
basis. 

^  Most  of  the  qualitative  results  obtained  herein  are  relative  to 

a  simulated  system  which  may  or  may  not  be  reasonable  relative  to  ship 
motions  problems.  An  experimental  or  analytical  investigation  of 
the  nature  of  cubic  frequency  response  funct ions  for  a  realistic  ship 
*  seakeeping  situation  is  recommended. 

The  "bottom  line"  in  all  seakeeping  predictions  is  usually  certain 
statistics  of  response  maxima,  and  only  when  this  is  possible  can  a 
^  complete  theory  of  seakeeping  be  claimed.  Efforts  should  thus  be  made 

to  develop  an  approach  to  the  probability  density  function  of  the 
maxima  of  a  cubic  system. 
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APPENDIX  A 

EXPECTED  VALUES  OF  PRODUCTS  OF 
ZERO  MEAN  GAUSSIAN  EXCITATION 

The  autocorrelation  function  or  second  moment  of  the  zero  mean 
Gaussian  excitation,  X(t),  is  defined  in  Equation  35  as: 

RxxW  =■  xTTT  m  -  t) 

By  definition  of  a  zero  mean  process: 

Xltl  =»  0  (A-l) 

From  the  development  of  Reference  25,  pages  82-85: 

X(t1)  X(t2)  •••  X(tn)  *  0  for  n  odd  (A-2) 


In  the  present  work  occasion  is  found  for  the  manipulation  of 
up  to  sixth  order  expectations.  Equations  35  of  the  text  and  Equa¬ 
tions  A-1  and  A-2  define  all  but  the  fourth  and  sixth  order  expecta¬ 
tions.  The  fourth  order  case  is  given  in  numerous  texts  including 
Reference  25,  and  may  be  written: 


X(t,)  X(t2)  XU,)  X(t4)  -  R^U,  -  t2)  Rxx(t3  -  tx) 

*  Rxx(t!  •  *3>  Rxx(t2  - 

+  RXX^£1  "  ^  RXX^2  "  C  3^ 


(A-3) 


where  the  zero  mean  assumption  has  been  utilized. 

The  sixth  order  expectation  as  derived  from  the  general  treat¬ 
ment  in  Reference  25  for  the  zero  mean  process  may  be  written  in 
present  notation  as: 
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x(tj)  xTt^7  x(t3)  xuj  x(t5)  x(tg) 

*  R  (t  -  t.)  R  (t,  -  t.  )  R  (t,  -  t,) 

XX  1  2  XX  3  4'  XX  5  6' 


+  RXX(tl  •  (2>  Rxx(t3  '  [5>  "xxK  ' 


+  R  (t  -  t.)  R  (t,  -  tc)  R  (t.  -  t,) 

XX  1  2  xx  3  6  XX  4  5 


*  Rxx(tl  -  *3>  Rxx(t2  -  '*>  Rxx(t5  '  r6> 


+  R  (t  -  tj  R  (t.  -  tc)  R  (t.  -  t,) 

XX  1  3  XX  2  5'  xx  4  6 


+  RXX^1  "  l3^  RXX^2  ■  t6^  Rxx^4  "  t5^ 
+RXX(tl  -  C4>  Rxx(t2  -  *3>  RXX(tS  -  »6> 
+  Rxx(tl  -  RXX(t2  •  *5>  Rxx(t3  -  V 

+  Rxx^l  ”  t4^  RXX^2  "  t()  Rxx^3  "  t5^ 

+  RXX^fl  "  t J  RXX^3  ’  c4^  Rxx^2  ‘  V 

+Rxx{tl  -  *5*  Rxx(t3  •  *2>  Rxx(t4  -  t6) 
+  RXX<tl  ■  ts>  Rxx^3  ■  *6*  Rxx^2  '  *4^ 

+  RXX(tl  "  V  RXX^3  "  th'>  Rxx^5  “  V 

+  RXX(tl  "  t6^  RXX^3  "  *5^  RXX^4  -  t2* 


+  RXX^1  ”  *6^  RXX^3  "  *2^  RXX^4  "  V 


(A- 3) 
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APPENDIX  B 

THE  n-DIMENS I ONAL  PARSEVAL  FORMULA 

The  n-dimensional  form  of  Parseval's  formula  finds  application 
in  the  present  work.  It  is  as  follows  (Reference  k) : 

VW  £n}  f2(tl't2*  •”  tn)dtldt2  •••  dtn 


(2*)" 


F1((i)1,u i2,  0Jn)  F2(aJl,w2,  •••  w^diOjdu^  •••  du)n 

(B-1) 


where  the  (*)  denotes  complex  conjugate  and  f j ( t j  •••)  and  Fj(tJ  *••) 
are  Fourier  Transform  pairs  defined: 

Fj  («*>!  ,u>2,  •••  un)  -  J  •••  |  fj(ti,t2,  tn)  Exp(-i  l  wrt(.)dt1dt2  •••  dtR 
fJ(W  **•  tn)  ***  \  Fj(W  —  wn)Exp(i  r^“rtr)dwl  ***  dwr 


(B-2) 
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A,  A. 


C. 

J 

»>> 

9,(1) 


«2<Vt2) 

93(tl,t2,t!) 


J 

9jk 

3jkt 
Gx  (u) 

G2  ,a)2^ 

G 3  (w^  »0J2  > u 3  ) 

j,k,i 

N 

P.P.r 

P  (j) 
yx  J 

P  (j ,k) 

yxx  J 

P  (j,M) 
yxxx  J ’  * 

R  (t) 

XX 

R  (t) 
yx 

Ryxx^tl’t2^ 
Ryxxx^ll ft2’ 


PRINCIPAL  NOTATION 

Constants 

Constants 

Constants 

Auxilliary  function  *Aa2+Ba+C 

n  n  n 

Linear  impulse  response 
Quadratic  impulse  response 
Cubic  impulse  response 
Oigital  impulse  response  (linear) 
Digital  impulse  response  (quadratic) 
Digital  impulse  response  (cubic) 
Linear  frequency  response  function 
Quadratic  frequency  response  function 
Cubic  frequency  response  function 
Time  increments 

Number  of  points  in  a  time  series 
Frequency  increments 
Second  order  periodogram 
Third  order  periodogram 
Fourth  order  periodogram 
Autocorrelation  function 
Cross  correlation  function 
Third  order  correlation 
Fourth  order  correlation 
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S  (w) 

XX 

Two  sided  spectrum,  excitation 

S  (oj) 

yy 

Two  sided  spectrum,  response 

S  (u) 
yx 

Cross  spectrum 

S  (nil  ,1*),) 
yxx  1  2 

Cross-bi -spectrum 

S  (<*),  .w-  ,u. ) 

yxxx  1  2  3 

Cross- tr i -spectrum 

'•'j 

Time 

V“’ 

One  sided  spectrum,  excitation 

u  (<*») 
yy 

One  sided  spectrum,  response 

w(k) 

FFT  of  squared  excitation 

x(t) 

Excitation 

4-» 

<  X 

Periodic  excitation 

X(k) 

FFT  of  excitation 

Y(t) 

Response 

Y(t) 

Response  to  periodic  excitation 

Y(k) 

FFT  of  response 

Z(k) 

FFT  of  cubed  excitation 

At 

Sampling  interval 

e.ej 

Phase  angles 

u,u. 

Circular  frequency 
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